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ABSTRACT 

We present a two-dimensional (2-D) fitting algorithm (Galfit, Version 3) with new capabilities to 
study the structural components of galaxies and other astronomical objects in digital images. Our 
technique improves on previous 2-D fitting algorithms by allowing for irregular, curved, logarithmic 
and power-law spirals, ring and truncated shapes in otherwise traditional parametric functions like 
the Sersic, Moffat, King, Ferrer, etc., profiles. One can mix and match these new shape features 
freely, with or without constraints, apply them to an arbitrary number of model components and 
of numerous profile types, so as to produce realistic-looking galaxy model images. Yet, despite the 
potential for extreme complexity, the meaning of the key parameters like the Sersic index, effective 
radius or luminosity remain intuitive and essentially unchanged. The new features have an interesting 
potential for use to quantify the degree of asymmetry of galaxies, to quantify low surface brightness 
tidal features beneath and beyond luminous galaxies, to allow more realistic decompositions of galaxy 
subcomponents in the presence of strong rings and spiral arms, and to enable ways to gauge the 
uncertainties when decomposing galaxy subcomponents. We illustrate these new features by way of 
several case studies that display various levels of complexity. 

Subject headings: galaxies: bulges — galaxies: fundamental parameters — galaxies: structure - 
techniques: image processing — techniques: photometric 



1. INTRODUCTION 

Images of astronomical objects store a wealth of infor- 
mation that encodes the physical conditions and fossil 
records of their evolution. Over the past decade, the 
ability of optical/near-infrared telescopes to resolve ob- 
jects improved by a factor of 10, and to detect faint 
surface brightnesses by at least 2 orders of magnitude. 
These advances now enable the study of highly intricate 
details on subarcsecond scales (e.g., nuclear star clus- 
ter, spiral structure, bars, inner ring, profile cusps, etc.), 
and extremely faint outer regions of galaxies. Moreover, 
new integral-field imaging capabilities blur the tradi- 
tional boundary of obtaining, analyzing, and interpreting 
imaging and spectroscopic data. Faced with the conver- 
gence in volume, quality, and multiwavelength datasets 
like never before, one of the main challenges toward mak- 
ing full use of the investments is developing sophisticated 
ways to extract information from the data to facilitate 
new science. 

1.1. Parametric and Non-Parametric Analysis 

Analyzing astronomical images is challenging because 
of the diversity in object sizes and shapes, and nowhere 
is it more difficult than for galaxies. Since the early era 
of photographic plates, one of the key methods for study- 
ing the light distribution of galaxies is to model it by us- 
ing analytic functions — a technique known as parametric 
fitting. This technique was first applied to galaxies by 
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Ide VaucouleursI (|1948f ) who showed that the light distri- 
bution of elliptical galaxies tended to follow a power-law 
form of exp (— r 1 / 4 ) . Subsequently, one of the break- 
throughs in our unde rstanding of gal axy structure and 
evolution came when iFreemar] ([1970ft showed that dy- 
namically "hot" stars in galaxies make up spheroidal 
bulges having a de Vaucouleurs light profile, whereas 
"cold" stellar components make up the more flattened, 
rotationally supported, exponential disk region. 

From that simple beginning, parametric fitting has 
been the mainstay for galaxy imaging studies, and ex- 
panding into many applications whenever the science 
calls for detailed and rigorous analysis. Among some of 
the examples, past investigations delved into the struc- 
tural parameters of disk galaxies (e.g.. Ide" Jong] [1996ft. 
the Tully-Fishe r relation (e.g iTullv k Fisherl Il977t 
iHinz et al.l 1200 3): iBedreeal et al.l [20061), the evolution of 
disky galaxies (ISimard et"aT1 120021 : iRavindranath et all 
12004 iBarden et alJl2005ft. the cosmi c evolution of galaxy 
morphology (e.g., iLillv et al.l [1998; Ma rleau fc Simardl 
119981 : lHathi et all 12009ft in both groundbased surveys 
and Hubble (Ultra-) Deep Fields (|Williams et al.l 119961: 
iBeckwith et al.ll2006ft . the morphological transformatio n 
of galaxies in cluster environments (e.g., [D ressier 1980ft. 
the fundamental plane of spheroids dDiorgovski &: Davisl 
119871: iDressler et al J IH387U Bender et al.l 1992ft. the red 



seque nce of galaxies (B ell et al J l2004bl fa nFtLer" et al.l 
l2007ft . morphological d issimilarities between s pheroidal 
galaxies and ellipticals ()Kormendvlll985l 11987ft . the cen 
tral structu r e of e a rly-type gal a xies 
Lauer et al.l 119951: iFaber et all ' 



d. 



(iKormendvl 119851 
Il997t iFerrarese et al 
2006bl lal: ILauer et al.ll2007ft an d implications for the for- 
mation of massive black hol es (|Ravindranath et al.ll2002l : 
IKormendv &: B ender 2009), black ho le vs. bulge re- 
lations ( |Kormendv fc Richstond I1995D a nd their evolu- 
tion (|Rix et all 120011: iPeng et al.l l2006aM bh. the "extra 
light" due to gas dissipation in galaxy centers ([Kormendvl 



fl99l iKormendv etaLl [20091: iHopkins et all l2008bllah 



quasar host galaxies (e.g., Hutchings et all 1984: 



McLeod fe Rie ko 1994; McLure et al.ll2000tlJahnke et all 



2004 ISanchez et all 120041: iKim et al.l l2008bD. gr avita- 
tional lensing of q uasar host galaxies ( Rix et all 120011 : 
iPeng et all l2006bT ). and t he clustering of dark matt er 
through weak lensing (e.g.. iHevmans etld"1l2006l , F2008). 

Since the original development of galaxy fitting nearly 
70 years ago, where the analysis was performed on 
one- dimensional (TP) surface brightness profile s (see 
alsolKqrmendvlll977l: lBursteidll97a lBorosodll98ll: iKent 
1985t lAndredakis fe Sanders! 11994 iMacArthur et al 



2001 . newer techniques have emerged t o directly analyze 



two-d im ensional (2-D) im ages (e.g. [Sh aw fe Gilmord 
1989t iMcLeod fc Riekd 11994 iBvun fc Freeman 



1995; 



1998; 



1999: 



1996; 



Ratnatu nga et aL. 
ISimard et all 120021: 



Moriondo et all 
1999t " 



19981: ISimard 



Wadadekar et ai- 
de Souza et all 12004 



Laurikainen et all 12004 iGadottil 120081 1 The bene 



fit of performing 2-D image analysis is to potentially 
make full use of all spatial information and to properly 
account for image smearing by the point-spread function 
(PSF). 

Even though 2-D analysis can be quite sophisticated, 
there are legitimate questions about whether it is more 
beneficial than 1-D for profile analysis. Proponents of the 

1- D technique are skeptical that perfect ellipsoid models 
are suitable to use for galaxies that show isophotal twists, 
or that are non-elliptical in shape. They note that not 
only is 1-D analysis more appealing because it is more 
straightforward to implement, the surface brightness pro- 
files serve as visual confirmation about the reality of fit- 
ting multiple component models. 

However, beneath the apparent simplicity there are a 
number of important subtleties to weigh. For instance, 
the decision about how to extract 1-D profiles is often 
not unique, nor are there strong reasons to prefer major 
or minor axis profiles, or a profile along some arc traced 
by spiral arms or isophote twists that result from the 
superposition of multiple components oriented at differ- 
ent angles. When symmetry is broken it is also unclear 
that there is an optimal or unique way to extract a 1-D 
profile, such as in irregular galaxies, overlapping galax- 
ies, and galaxies with double nuclei. Another factor to 
consider is that the process of extracting 1-D profiles re- 
duces spatial information content: in many situations, 
a bulge, disk, and bar can all have different axis ratios, 
position angles (PAs), and profiles that help to break 
model degeneracies, but this information is lost when the 
data are collapsed into 1-D. Lastly, for compact galax- 
ies, 1-D profile fitting cannot properly correct for image 
smearing by the PSF because 1-D profile convolution is 
not mathematically equivalent to convolution in a 2-D 
image. While some of the above concerns also affect 2- 
D analysis (i.e. irregular galaxies), most others benefit 
from treatment using 2-D techniques. When it comes to 
judging which models are more plausible, there are few 
diagnostics more discerning than a moment's glance at 

2- D models and residual images; a good fit in 2-D always 
means that 1-D profiles are necessarily a good fit. Propo- 
nents of 2-D analysis therefore believe that the benefits 
outweigh the drawbacks. Moreover, many drawbacks can 
be mitigated by breaking free from axisymmetry in 2-D 
analysis, which is the purpose of this study to show. 



In the box of tools for morphology analysis, a compli- 
mentary approach is non-parametric analysis. While we 
do not use non-parametric methods in this study, it is 
useful to understand the conceptual differences between 
the two approaches. We thus provide a brief overview. 
In contrast to function fitting, the non-parametric ap- 
proach does not involve deciding what functional form to 
use or how many. One method is to de compose an im- 
age into "shapelets " or "wavelets" (e.g., Rcfrcgicrl 120031 : 
iMassev et all 12004). which is analogous to taking a 2- 
D Fourier transform of an image using mathematically 
orthogonal basis functions. The main conceptual differ- 
ence with parametric fitting is that the shapelet basis 
functions do not represent physical subcomponents of a 
galaxy. Moreover, the power spectrum of the basis func- 
tions is quite useful for diagnosing the degree of galaxy 
distortions. There are also other non-parametric tech- 
niques (e.g., lAbraham et allll994 iRudnick fc Rb3 ll998t 
IConselice et al.ll2000t lLotz et alll2004D . To measure con- 
centration non-parametrically, one way is to compare 
fluxes within apertures of different radii; whereas to mea- 
sure asymmetry one can rotate an image by 180° and 
subtract it from the original image and measure the resid - 
uals (e.g., lAbraham et all 11994 IConselice et all [2 000). 
Towar d t he same goals, t wo studies. lAbraham et all 
( 2003) and lLotz et al.l (|2004D , introduce the Gini index to 
measure the concentration of a galaxy by comparing the 
relati ve distribution of pixel flux values within a certain 
area. lLotz et ail (|2004D also introduce a method for mea- 
suring asymmetry through the M20 parameter, which is 
the second-order moment of the brightest 20% of the a 
galaxy's flux. 

The application of non-parametric ana lysis has mostly 
been to quantify gala xy mergers (e.g., IConselice et all 
l2003t lLotz et al.ll2008l ). These techniques are generally 
much simpler to implement than parametric fitting and 
have a strong virtue that no assumptions are made about 
the galaxy profiles and shapes. The tradeoff is that the 
techniques often do not deal with image smearing by the 
PSF and different sensitivity thresholds between different 
surveys. Consequently, one has to take particular care 
to compare compact with extended objects, measured in 
different apertures, or me asured from i mages of different 
surface brightness depths (Liskcr 2001) . One also should 
guard against contamination by intervening galaxies or 
stars because the techniques do not have a rigorous way 
to separate overlapping objects. For separating objects, 
extracting structural components of a galaxy, and ex- 
trapolating galaxy wings well into the background noise, 
there are few, if any, alternatives to parametric analysis 
that are more rigorous. 

When comparing the merits of non-parametric and 
parametric analysis, the idea of ellipsoid models in para- 
metric analysis is sometimes considered to be a weakness, 
because galaxies, after all, are not perfect ellipses in pro- 
jection. However, it is worth pointing out that the no- 
tion of there being a global average size inherently implies 
comparison against some kind of approximate shape. In- 
deed, even in non-parametric techniques, to measure a 
size in an 2-D image, one assumes a basic shape either 
explicitly (through using aperture photometry) or implic- 
itly (through calculating flux moments, which requires a 
center to be defined). An ellipsoid is one of the sim- 
plest and most natural low-order shapes against which 
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all galaxies can be compared, especially for measuring 
an average size. This notion is useful: deviations from 
the basic ellipsoid shape can then be considered as higher 
order modifications, even for highly irregular galaxies. 

Nevertheless, there are many situations where it is de- 
sirable to use models that deviate from ellipsoid shapes. 
Contrary to the common practice, there are numerous 
ways to break from axisymmetry. However, the harder 
challenge is to devise a scheme that is intuitive to grasp 
and well motivated. Breaking free from axisymmetry al- 
lows for other interesting science applications, including 
a promising new way to quantify asymmetry. 

1.2. The Next Generation of Parametric Imaging 
Fitting 

In this study we present, as a proof of concept, new ca- 
pabilities in 2-D image fitting that progress beyond the 
limitations of traditional parametric fitting models. One 
key aspect of our approach is to first identify a mini- 
mum basis set of features that spans the range of galaxy 
morphologies and shapes. From experience, we deter- 
mine those four "basic elements" to be bending, Fourier, 
coordinate rotation, and truncation modes. Secondly, 
one of the main reasons why parametric fitting is use- 
ful is that the profile parameters are intuitive to grasp 
(e.g., concentration index, effective radius, total lumi- 
nosity, etc.). Therefore, another key requirement is that 
the traditional profile parameters must retain their origi- 
nal, intuitive, meaning even under detailed shape refine- 
ments, and even under such extreme cases as irregular 
galaxies. This can be accomplished if the basic premise 
starts with the traditional ellipsoid function, on top of 
which one can add perturbations, rotations, irregulari- 
ties, and curvature. This is possible because of the fact 
that simple ellipsoid fits are a reasonable way to quan- 
tify global average properties, and other details can be 
considered to be higher order perturbations that may be 
of other practical interest. 

As we attempt to demonstrate, combining just the four 
basic mophology elements can quickly yield a dizzying ar- 
ray of possibilities for fitting galaxies. The end result can 
look highly "realistic." Indeed, it is now possible to fit 
many spiral galaxies, asymmetric tidal features, irregu- 
lar galaxies, ring galaxies, dust lanes, truncated galaxies, 
arcs, among others (though, certainly, there are limita- 
tions). However, it is important to realize that "being 
possible" often does not mean "being necessary" or "be- 
ing practical." Necessity ought to be judged in the sci- 
entific context of whether it is worth the extra effort to 
obtain diminishing returns. For instance, to measure to- 
tal galaxy luminosity, it is often unnecessary to fit high 
order Fourier modes or spiral rotations. For many sci- 
ence studies interested in global parameterization, often 
a single ellipsoid component would suffice. It is there- 
fore important to always let the science determine what 
kind of analysis is required, rather than to use the new 
capabilities in the absence of a clearly defined goal. Hav- 
ing provided some foregoing disclaimers, some of the key 
scientific reasons motivating the new capabilities are to: 

• Quantify global asymmetry or substructure asym- 
metry. 

• Quantify bending modes for weak-lensing applica- 



tions, or fit arcs in the image plane for strong grav- 
itational lenses. 

• Obtain more accurate substructure decomposition 
in the presence of bars, spirals, rings, etc. 

• Obtain more accurate global photometry. 

• Quantify profile deviation in inner or outer regions 
of a galaxy, such as disk truncations, deviations 
from a Sersic function, etc. 

• Extract parametric information to the limits im- 
posed by resolution, signal-to-noise, and other 
small scale fluctuations. 

• Quantify model-dependent errors in the decompo- 
sition. 

We thus begin by giving an overview of the Galfit 
software (§|2]). Then, we introduce the radial profile func- 
tions that one can use (§0, and illustrate how symmet- 
ric and asymmetric shapes can be generated by modify- 
ing the coordinate system in (§ |4j) . Next, we introduce 
a new capability that allows for radial profile trunca- 
tion (§ [5]). Enabling all the capabilities may result in 
extremely complex galaxy shapes, the interpretation of 
which may give concerns to those new to the analysis. 
Therefore, we discuss the interpretations and model de- 
generacies of the parameters in § [6] We then apply these 
new features to real galaxies in § [7J followed by conclud- 
ing remarks (§[8|). 

2. THE 2-D FITTING PROGRAM Galfit 

This st udy builds on an existing algorithm named 
Galfit 5 (jPeng et all [2002) , which is a 2-D paramet- 
ric galaxy fitting algorithm, in the same spirit as 
other wi dely used 2-D image-fitting algor ithms (e.g., 
( GIM2D: IShriardl [1998L f Simard et all I2002T) : (BUDDA: 
Ide Souza et al.ll2004( )). Galfit is a stand-alone program 
written in the C language, and can be run on most mod- 
ern operating systems. To read and produce FI TS im- 
ages, Galfit calls upon the CFITSIO package (Pence 
1999). Galfit is designed to allow for complex im- 
age decomposition tasks: by allowing for an arbitrary 
number and mix of parametric functions (Sersic, Moffat, 
Gaussian, exponential, Nuker, etc.), it can simultane- 
ously fit any number of galaxies and their substructures. 
It is possible to use Galfit for both interactive anal- 
ysis and galaxy surveys where complete automation is 
required. However, automation requires the use of an 
external "wrapping" algorithm written by the user that 
takes care of both the pre-processing (object identifica- 
tion, initial parameter estimation) and post-processing 
(extracting and tabulating fitting parameters) of the fit- 
ting results. 

2.1. xt an d Error Analysis 

Galfit is a non-linear least-squares fitting algorithm 
that uses the Levenberg-Marquardt technique to find the 
optimum solution to a fit. The Levenberg-Marquardt 
algorithm is currently the most efficient one for searching 
large parameter spaces, allowing for the possibility to fit 

5 http: / /users. obs.carnegiescience.edu/peng/work/galfit/galfit. html 
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complex images with multiple components and a large 
number of parameters. Galfit determines the goodness 
of fit by calculating \ 2 and computing how to adjust the 
parameters for the next step. It continues to iterate until 
the x 2 n ° longer decreases appreciably. The indicator of 
goodness of fit is the normalized or reduced \ 2 > xt '■ 

2 1 (/data(z, y) - /modcl(z, V)) 2 

where 

m 

/model (a;, y) = y;ai •••"«)• ( 2 ) 

A^jof is the number of degrees of freedom in the fit; nx 
and ny are the x and y image dimensions; and /data(£, y) 
is the image flux at pixel (x,y). The f mo dei(x, y) is the 
sum of m functions of f„(x, y; a\...a n ), with n free pa- 
rameters (ai...a n ) in the 2-D model. The uncertainty as 
a function of pixel position, a(x,y), is the Poisson error 
at each pixel, which can be provided as an input image. 
If no cr-image is given, one is generated based on the 
gain and read-noise parameters contained in the image 
header. Pixels in the image marked as being bad do not 
enter into the calculation of \ 2 ■ 

In the Levenberg-Marquardt algorithm, the minimiza- 
tion process involves computing a Hessian matrix, which 
is closely relat ed to the covarian ce matrix of the param- 
eters (e.g., see iPress et al.l[l992l ). The covariance matrix 
is then directly related to the formal uncertainty in the 
fitting parameters that Galfit reports. However, the 
usefulness of the formal uncertainty is limited to ideal 
situations where the fluctuations in the residual image 
are only due to Poisson noise after removing the model. 
This situation is mostly realized in idealized situations, 
such as image simulations. In real images, the residuals 
are due to structures like stars and galaxies that are not 
fitted, flat-fielding errors, and imperfect functional match 
to the data. These factors cause formal uncertainties re- 
ported in numerical fits to be only lower estimates. In 
image fitting, more realistic uncertainties are necessarily 
obtained by other processes, such as comparing fit results 
based on different assumptions about the model rather 
than through a formal covariance matrix. 

In summary, the three images Galfit takes as input 
to calculate least squares are the data, a cr-image, and an 
optional bad pixel mask. To account for image smearing 
by the PSF, Galfit will also require a PSF image. 

2.2. Accounting for Telescope Optics and Atmospheric 

Seeing 

The wavefront of light from distant sources is always 
perturbed by the act of producing an image, distortions 
due to imperfect optics, and sometimes by the Earth's 
atmosphere, resulting in some blurring. To accurately 
compare the intrinsic shape of an object with a model, 
image blurring must be taken into account. In image fit- 
ting this is often done by convolving a model image with 
the input PSF before comparing with the data. The pro- 
cess of performing convolution is mathematically rigor- 
ous, but the actual implementation has several subtleties. 



One consideration is the computation speed, as the 
process of convolving a model is frequently the most time 
consuming part of parametric fitting. The trade-off is 
that the smaller the convolution region the faster the 
computation time, but also the less accurate. To achieve 
a compromise, Galfit allows the user to decide on the 
size of the convolution region. This gives the flexibility 
for one to hone in on a solution quickly before trying to 
obtain higher accuracy in the final step. 

Another important issue to consider is whether to con- 
volve each component separately or all of them together 
just once in the final image. This is an important con- 
sideration because even though the model functions are 
analytic, they are resampled by discrete pixel grids, re- 
sulting in a "pixellated" profile instead of one that is 
infinitely smooth. If an intrinsic model is sufficiently 
sharp, the curvature may not be critically subsampled 
by the pixels prior to convolution, regardless of whether 
the recorded data are Nyquist sampled. The resulting 
profile after image convolution therefore can depend very 
sensitively on how the model is centered on a pixel. If 
such a model is created off-center, pixellation effectively 
broadens out the model ever so slightly more than normal 
once convolution is applied, but the effect is noticeable in 
high-contrast imaging studies. Therefore, the better way 
to deal with "pixellation broadening" is to convolve each 
model component individually rather than the entire im- 
age at once. To do so, Galfit creates every model on a 
pixel center; the pixel fluxes near the center of the mod- 
els are integrated over the pixel area adaptively. Then 
to effect an off-centered model, Galfit makes use of 
the convolution theorem by shifting the PSF by the re- 
quired amount before convolving it with the model. This 
process circumvents the problem of artificial pixellation 
broadening because whereas the model core region may 
not be sufficiently resolved, the PSF ought to be 6 . 

Shifting the PSF, however, can be quite problematic 
when it is marginally Nyquist sampled, or if the diffrac- 
tion patterns are not critically sampled. Accurate shift- 
ing of the PSF is of basic importance in high contrast 
imaging studies. For instance, in the case of studying 
active galaxies with a strong central point source, issues 
of contrast, resolution, and sampling all conspire to make 
the PSF fitting crucial to deriving a reliable host model. 
In such situations, the standard interpolation techniques 
(e.g., linear or spline) tend to broaden out the PSF core, 
so they are only accurate in the extended outskirts where 
the gradient is shallow. One alternative is to interpolate 
using the sine kernel, which is theoretically the perfect 
interpolation kernel for critically sampled images, and 
preserves the intrinsic width of the data. However, signif- 
icant "ringing" appears around sharp features (i.e. PSF 
or galaxy core). This effect can be nearly as bad on a 
fit as pixellation broadening. An improved solution is 
to taper the wing of the sine kernel using a windowing 
function (e.g., Lanczos), but the ringing often may still 
be quite large beyond the PSF core, which must be fur- 
ther suppressed. 

Galfit seeks a compromise by using a hybrid scheme 
where the interpolation in the PSF core is done by us- 
ing a sine kernel with a Kaiser window function so as 

6 If the PSF is not resolved then the convolution process will not 
be accurate regardless of the technique. 
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# IMAGE and GALFIT 

A) gal. fits 

B) imgblock.f its 

C) none 

D) psf.fits # 

E) 1 

F) none 

G) none 

H) 1 200 1 

I) 100 100 
J) 20.000 
K) 1.000 1.000 
Q) both 
P) 



CONTROL PARAMETERS 

# Input data image (FITS file) 

# Output data image block 

# Sigma image name (made from data if blank or "none") 

# Input PSF image and (optional) diffusion kernel 

# PSF fine sampling factor relative to data 

# Bad pixel mask (FITS image or ASCII coord list) 

# File with parameter constraints (ASCII file) 
100 # Image region to fit (xmin xmax ymin ymax) 

# Size of the convolution box (x y) 

# Magnitude photometric zeropoint 

# Plate scale (dx dy) [arcsec per pixel] 

# Display type (regular, curses, both) 

# Choose: 0=optimize, l=model, 2=imgblock, 3=subcomps 



INITIAL FITTING PARAMETERS 

For component type, the allowed functions are: 

sersic, expdisk, edgedisk, devauc, king, nuker, psf, 
gaussian, moffat, ferrer, and sky. 

Hidden parameters will only appear when they're specified: 
Bn (n=integer, Bending Modes). 
CO (diskiness/boxiness) , 
Fn (n=integer, Azimuthal Fourier Modes) . 

R0-R10 (coordinate rotation, for creating spiral structures). 
To, Ti, T0-T10 (truncation function). 



par) par value(s) fit toggle(s) # parameter description 



# Component number : 1 



0) 
1) 
3) 
4) 
5) 
9) 
10) 
Ti) 



sersic3 / 
50.0000 50.0000 



15.0000 
30.0000 
4.0000 
0.7000 
-30.0000 
2 



# Component type 
1 1 # Position x, y 

# Surface brghtnss outer R_break [mag/arcsec~2] 

# R_e (effective radius) [pix] 

# Sersic index n (de Vaucouleurs n=4) 

# Axis ratio (b/a) 

# Position angle (PA) [deg: Up=0, Left=90] 

# Inner truncation by component number(s) 



F5) 0.1500 20.0000 1 1 # Azim. Fourier mode 5, amplitude, & phase angle 



# Component number : 2 
TO) radial 



TI) 45.0000 
T4) 8.0000 
T5) 5.0000 
T9) 0.7000 
T10) 45.0000 



# Truncation type (radial, length, height) 

45.0000 1 1 # Position x, y 

1 # Break radius (99'/, normal flux) [pixels] 

1 # Softening length (1% normal flux) [pixels] 

1 # Axis ratio (optional) 



# Position angle (optional) [deg: Up=0, Left=90] 
Fl) 0.6000 20.0000 1 1 # Azim. Fourier mode 1, amplitude, & phase angle 
B2) -5.000e+00 1 # Bending mode 2 amplitude 



# Component number : 3 
0) 
1) 
3) 
4) 
5) 



sersic 

150.0000 50.0000 
7.0000 



15.0000 
2.0000 
9) 0.5000 
10) 0.0000 



# Component type 
1 1 # Position x, y 

# Integrated magnitude 

# R_e (effective radius) [pix] 

# Sersic index n (de Vaucouleurs n=4) 

# Axis ratio (b/a) 

# Position angle (PA) [deg: Up=0, Left=90] 



R0) power 
Rl) 0.0000 
15.0000 
180.0000 



R2) 
R3) 



R4) 0.3000 
R9) 10.0000 
R10) 45.0000 
Fl) 0.3000 
F5) 0.1000 



1 
1 
1 
1 
1 
1 

45.0000 
90.0000 



# PA rotation f unc . (power, log, none) 

# Spiral inner radius [pixels] 

# Spiral outer radius [pixels] 

# Cumul . rotation out to outer radius [degrees] 

# Asymptotic spiral powerlaw 

# Inclination to L.o.S. [degrees] 

# Sky position angle 

1 # Azim. Fourier mode 1, amplitude, & phase angle 

1 # Azim. Fourier mode 5, amplitude, & phase angle 



Fig. 1. — Example of an input file. The object list is dynamic and can be extended as needed. Each model is modified by a mix of higher 
order Fourier modes, bending modes, truncation, or spiral structure. These parameters produce the models shown in Figure[2] 
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Fig. 2. — Shapes produced by parameters in the Galfit input 
file of Figure [T] Left: a Sersic light profile modified by a single 
Fourier mode m = 5, creating the star shape. It is truncated in 
the inner region by a truncation function, which is modified by a 
bending mode m = 2, with a lopsided Fourier mode of m = 1. 
Right, a Sersic light profile with Fourier modes m = I and m = 5 
is modified by a coordinate rotation function to create a lopsided, 
multi-armed, spiral structure. 



to faithfully preserve the width, but a bicubic spline 
interpolation is used in the wings. The result of this 
scheme is that for a Gaussian having a full width at half- 
maximum (FWHM) of 2 pixels, the interpolation is ac- 
curate to 0.1% in the center, and 0.03% at the distance 
of the FWHM relative to the peak (or 1% relative to 
the local flux). For oversampled PSFs, the interpola- 
tion is even more accurate. Compared to bicubic spline 
interpolation, our scheme is about 20 times more accu- 
rate. From a more practical standpoint, the mismatch 
in the PSFs between data taken using the Hubble Space 
Telescope (HST) imaging cameras and synchronously ob- 
served PSFs is rarely better than 3% in the core. For 
all practical purposes, our interpolation scheme there- 
fore will more than suffice for the most demanding high 
contrast studies of quasar host galaxies at high redshift. 

When the data are undersampled, convolution of the 
model can still be done correctly if the convolution PSF 
provided to Galfit is either critically sampled or over- 
sampled. In this situation, Galfit will generate a model 
on a finer grid, convolve it with the PSF, then bin the 
result down to the resolution of the data for comparison. 
One way for users to obtain an oversampled PSF com- 
pared to the data is to dither the PSF observations by 
fractional pixels. Another way is to numerically recon- 
struct a more oversampled PSF star by extracting multi- 
ple stars from the data image itself (e.g., via DAOphot, 
IStetsonlll987t ). 

However, lastly, we note that when the data and the 
convolution PSF are both undersampled (i.e. with PSF 
FWHM < 2 pixels), convolution cannot be done accu- 
rately. In such a situation, for the purpose of image 
fitting, it is often better to broaden out the data and the 
PSF to critical samplin g than to perform the analysis in 
the original resolution (|Kim et al.ll2008al ). 

2.3. The Concept of a Model Component 

Using the new features, each model can take on a 
shape that is completely unrecognizable from a tradi- 
tional ellipsoid shape. It is therefore necessary to clarify 
what constitutes a single model component. In Gal- 
fit, each model component is referred to by the name 
of the surface brightness profile, just as it is standard 



practice to call something a Sersic, Gaussian or expo- 
nential component in traditional models. As implied by 
this notion, no matter how complex the shape, the flux 
declines monotonically (unless modified by a truncation 
function, Section [5| from a peak in every radial direc- 
tion in a non-rotating frame, or along an arc in a rotat- 
ing frame, strictly following the functional form specified 
by the user. The radial profile parameters are mathe- 
matically decoupled from the azimuthal shape because 
the radial profile functions are self-similar in the expres- 
sion of the radius parameter, i.e. with powers of (r/r e ), 
whereas the complex azimuthal shapes are obtained by 
simply stretching the coordinate metric into more exotic 
grids than the standard Cartesian grid. This idea is in 
fact implicit in all 2-D image-fitting algorithms, where 
the axis ratio parameter, q, turns a circular profile into 
an ellipse by compressing the coordinate axis along one 
direction, even though the functional form of the profile 
remains the same in every direction. In the same man- 
ner, the definition of a scale or effective radius in a com- 
ponent, no matter how complex the shape, corresponds 
closely to that of the best-fitting ellipse in the direction 
of the semi-major axis. 

Figure[T]demonstrates how Sersic profiles can be modi- 
fied by bending modes, Fourier modes, and a spiral rota- 
tion function in Galfit — the results of which are shown 
in Figure O In the example, there are only two Sersic 
model components, despite the appearance of numerous 
parameters: both the "radial" and "power" functions are 
modifications to the Sersic profiles. Furthermore, the 
Fourier and bending modes can modify the Sersic pro- 
files, or modify the modifiers to the light profiles. Each 
radial surface brightness profile has a single peak and the 
flux decline is monotonic radially (unless truncated by a 
truncation function called "radial" in Figure [I} or in a 
rotating coordinate system (called "power" in Figure [IJ . 
Therefore, for each component, it is still meaningful to 
talk about, for example, an "average" light profile (e.g., 
Sersic), with an average Sersic concentration index n — 
no matter what the galaxy may look like azimuthally. In 
this manner even irregular galaxies can be parameterized 
in terms of their average light profile. When the average 
peak of an irregular galaxy is not located at the geomet- 
ric center, it has a high-amplitude m — 1 Fourier mode 
(i.e. lopsidedness) . 

In such a way, no matter how complex the azimuthal 
shape, interpreting the surface brightness profile parame- 
ters is just as straightforward as the traditional ellipsoid. 

3. THE RADIAL PROFILE FUNCTIONS 

The radial profile functions describe the intensity fall- 
off of a model away from the peak, such as the Sersic, 
Nuker, or exponential models, among others. For exam- 
ple, early-type galaxies typically have steep radial pro- 
files whereas late-type galaxies have shallower intensity 
slope near the center. The rate of decline is governed by 
a scale-length parameter. The radial profile is often of 
primary interest in galaxy studies from the standpoint 
of classification, and because the exact functional form 
may have some bearing on the path of galaxy evolution. 
In Galfit the radial profile can have the following func- 
tional forms, which are some of the most frequently seen 
in literature. 
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Fig. 3. — The Sersic profile, where r e and S e are held fixed. 
Notice that the larger the Sersic index value n, the steeper the 
central core, and more extended the outer wing. A low n has 
a flatter core and a more sharply truncated wing. Large Sersic 
index components are very sensitive to uncertainties in the sky 
background level determination because of the extended wings. 



The Sersic Profile The Sersic power law is one of 
the most frequently used to study galaxy morphology, 
and has the following functional form: 



S(r) = E e exp 




(3) 



S e is the pixel surface brightness at the effective radius 
r e . The parameter n is often referred to as the concen- 
tration parameter. When n is large, it has a steep inner 
profile and a highly extended outer wing. Inversely, when 
n is small, it has a shallow inner profile and a steep trun- 
cation at large radius. The parameter r e is known as the 
effective radius such that half of the total flux is within 
r e . To make this definition true, the dependent variable 
k is coupled to n; thus, it is not a free parameter. The 
classic de Vaucouleurs profile that describes a number 
of galaxy bulges is a special case of the Sersic profile 
when n = 4 (corresponding to ft = 7.67). As explained 
below, both the exponential and Gaussian functions are 
also special cases of the Sersic function when n = 1 and 
n = 0.5, respectively. As such the Sersic profile is a 
common favorite when fitting a single component. 
The flux integrated out to r = oo for a Sersic profile 

is: 

F tat = 2TTr 2 e Z e e K n K - 2n r(2n)q/R(C a ; m). (4) 

The term i?(Co; mi) is a geometric correction factor when 
the azimuthal shape deviates from a perfect ellipse. As 
the concept of azimuthal shapes will be discussed in Sec- 
tion 21 we will only comment here that R(Cq; mi) is sim- 
ply the ratio of the area between a perfect ellipse with 



the area of the more general shape, having the same 
axis ratio q and unit radius. The shape can be modi- 
fied by Fourier modes (m being the mode number) or 
diskiness/boxiness. For instance, when the shape is mod- 
ified by diskiness/boxiness, R(Cq) has an analytic solu- 
tion given by: 



R(Co) 



HCo + 2) 



4/3(l/(Co + 2),l + l/(C + 2))' 



(5) 



where /? is the Beta function. In general, when the 
Fourier modes are used to modify an ellipsoid shape, 
there is no analytic solution for R(mi), and so the area 
ratio must be integrated numerically. 

In Galfit, the flux parameter that one can use for the 
Sersic function is either the integrated magnitude mtot 
or some kind of surface brightness magnitude, for exam- 
ple at the center (/io), at the effective radius (/x e ), or at 
the break radius (//break) for truncated profiles (see Sec- 
tion [5]). The integrated magnitude follows the standard 
definition: 



m to t 



-2.51o gl 



exp 



mag zpt, 



(6) 



where t cxp is the exposure time from the image header. 
Each Sersic function can thus potentially have 7 classical 
free parameters in the fit: xq, yo, mtot, r e , n, q, and 
Qpa- The non-classical parameters, Co, Fourier modes, 
bending modes, and coordinate rotation may be added as 
needed. There is no restriction on the number of Fourier 
modes, and bending modes, but each Sersic component 
can only have a single set of Cq and coordinate rotation 
parameters (see Section [4] for details). 

The Exponential Disk Profile The exponential 
profile has some historical significance, so Galfit is ex- 
plicit about calling this profile an exponential disk, even 
though an object that has an exponential profile need not 
be a classical disk. Historically, an exponential disk has 
a scale length r s , which is not to be confused with the ef- 
fective radius r e used in the Sersic profile. For situations 
where one is not trying to fit a classical disk it would be 
less confusing nomenclature-wise to use the Sersic func- 
tion with n = 1, and quote the effective radius r e . But 
because the exponential disk profile is a special case of 
the Sersic function for n = 1 (see Figure [3|) , there is a 
relationship between r e and r s , given by 

r e = 1.678r s (For n=l only). (7) 

The functional form of the exponential profile is 



£(r) = So exp 
and the total flux is given by 

F to t -27rr s 2 S g/i?(C ;m). 



(8) 



(9) 



The 6 free parameters of the profile are: xq, yo, mtot, r s , 
6>pa, and q. 

The Gaussian Profile The Gaussian profile is an- 
other special case of the Sersic function with n = 0.5 (see 
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Fig. 4. — The modified Ferrer profile. The black reference curve 
has parameters r out = 100, a = 0.5, /3 = 2, and So = 1000. The 
red curves differ from the reference only in the a parameter, as 
indicated by the red numbers. Likewise, the green curves differ 
from the reference only in the j3 parameter, as indicated by the 
green numbers. 
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Fig. 5. — The empirical King profile. The black reference curve 
has parameters r c = 50, rt = 100, a = 2, and So = 1000. The 
red curves differ from the reference curve only in the a parameter, 
as indicated by the red numbers. Likewise, the green curves differ 
from the reference only in the r c parameter, as indicated by the 
green numbers. 



FigureEJ, but here the size parameter is the FWHM in- 
stead of r e . The functional form is 

£(r) = £ expfe), (10) 
and the total flux is given by 

F to t = 2w 2 E g/i?(C ;m), (11) 

where FWHM = 2.354<t. The 6 free parameters of the 
profile are: Xq, yo, m-tot: FWHM, q, and #pa- 

The Modifie d Ferrer Profil e The Ferrer profile 

(Figure [4j iBinnev fc Tremaind Il987f) has a nearly flat 
core and an outer truncation. The sharpness of the trun- 
cation is governed by the parameter a, whereas the cen- 
tral slope is controlled by the parameter /?. Because of 
the flat core and sharp truncation behavior, historically 
it is often used to fit galaxy bars and "lenses." The pro- 
file 

S(r) = S„(l-(r/r out )^)" (12) 

is only defined within r < r out , beyond which the func- 
tion has a value of 0. The 8 free parameters of the Ferrer 
profile are: xq, yo, central surface brightness, r out , a, /3, 
q, and 6pa- 

It is worth mentioning that a Sersic profile with low 
index n < 0.5 has similar profile shapes, thus it is often 
used instead of the Ferrer function. 
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Fig. 6. — The Moffat profile. The black reference curve has 
parameters n = 2, FWHM = 20, and S = 1000. The other 
colored lines differ only in the concentration index n, as shown by 
the numbers. The dashed line shows a Gaussian profile of the same 
FWHM. 



The Empirical (Modified) King Profile The em- 
pirical King profile (Figure [5J is often used to fit the 



light profile of gl obular clusters. It has the following 
form (|El^[l99l : 
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Fig. 7. — The Nuker profile. The black reference curve has parameters r;, = 10, a = 2, /3 = 2, 7 = 0, and If, = 100. For the other colored 
lines, only one value differs from the reference, as shown in the legend. 
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(1 + (r t /r c )2)V« 
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(l + Cr/re) 2 ) 1 /- (1+ (r t /r c ) 2 )V° 



.(13) 



The standard empirical King profile has a power law with 
index a = 2. In Galfit, a can be a free parameter. 
In this model, the flux parameter to fit is the central 
surface brightness, /^o, expressed in mag arcsec -2 (see 
Equation |20|) . The other free parameters are the core 
radius (r c ) and the truncation radius (r t ), in addition 
to the geometrical parameters. Outside the truncation 
radius, the function is set to 0. Thus, the total number 
of classical free parameters is 8: xq, yo, fj,Q, r c , rt, a, q, 
and 6pA- 

The Moffat Profile The profile of the HST WFPC2 
PSF is well described by the Moffat func tion (Figure E3) . 
Other than that, the Moffat function (jMoffatl fl969l) is 
less frequently used than the above functions for galaxy 
fitting. The functional profile is 



E(r) 



and the total flux is given by 

(n - l)R(C ;m)' 



-ft.nt — 



(14) 



(15) 



In Galfit the size parameter to fit is the FWHM, where 
the relation between and FWHM is 



FWHM 

2V2V« - l' 



(16) 



The 7 free parameters are: xq, yo, m to t (i-e. total mag- 
nitude, instead of fio) FWHM (instead of r^), the con- 
centration index n, q, and 9pa- 

The Nuker Profi le Th e Nuke r profile (Figure [7]) was 
introduced by lLauer et al.l ()1995[ ) to fit the central light 
distribution of nearby galaxies, and it has the following 
form: 



J(r) 



h 2 







1 



(17) 



The flux parameter to fit is /if,, the surface brightness of 
the profile at n,, which is defined as 



A*fc 



-2.5 lo gl 



h 



t cx _ p AxAyJ 



mag zpt , 



(18) 



where t cxp is the exposure time from the image header, 
and Ax and Ay are the platescale in arcsec. The Nuker 
profile is a double power law, where (in Equation 1171) /3 
is the outer power law slope, 7 is the inner slope, and a 
controls the sharpness of the transition. The motivation 
for using this profile is that the n uclei of many gala xies 
appear to be fit well in 1-D (see lLauer et al.(ll995T ) by 
a double power law. However, caution should be exer- 
cised when using this function because, for example, a 
low value of a (a < 2) can be mimicked by a combina- 
tion of high 7 and low j3 (compare Figure [TJc with the 
other two panels), which presents a serious potential for 
degeneracy. In all there are there are 9 free parameters: 
xo, ?/o, Mb, n, a, (3, 7, q, and 6> PA - 

The Edge-On Disk Profile Both the Sersic (Equa- 
tion [3J and exponential disk profile (Equation [8]) are 
merely empirical descriptors of a galaxy light profile. 
However, for edge-on disk galaxies, there is a more 
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physically motivated light profile: under the assump- 
tion that the disk component is locally isothermal and 
self-gravitating, the light pro file distribution is given by 
(jvan der Kruit fc Searlelll981l ): 



E(r,/i) = S 



(19) 



where So is the pixel central surface brightness, r s is 
the major-axis disk scale length, h s is the perpendicular 
disk scale height, and K\ is a Bessel function. The flux 
parameter being fitted in Galfit is the central surface 
brightness: 



Ai = -2.5 log 



10 



cxp 



§^) + mag Z pt. (20) 



Note that if the disk is oriented horizontally the coor- 
dinate r is the x-distance (as opposed to the radius) of a 
pixel from the origin. There are 6 free parameters in the 
profile model: x , yo, Mo> r s , h s , and P a- 

The PSF Profile For unresolved sources, one can 
fit pure stellar PSFs to an image (as opposed to func- 
tions with narrow FWHM convolved with the PSF). The 
PSF function is simply the convolution PSF image that 
the user provides, so there is no prescribed analytical 
functional form. This is also the only profile that is not 
convolved in Galfit. The PSF has only 3 free param- 
eters: xq, yo, and TOtot- Because there is no analytical 
form, the total magnitude is determined by integrating 
over the PSF image and assuming that it contains 100% 
of the light. If the PSF wing is vignetted, there will be 
a systematic offset between the flux Galfit reports and 
the actual value. 

If one wants to fit this "function," it is important to 
make sure that the input PSF is close to, or super-, 
Nyquist sampled. The PSF interpolation used in shifting 
is done by a sine function with a Kaiser window, which 
can preserve the widths of the PSF even under subpixel 
shifting. This is in principle better then spline interpo- 
lation or other high-order interpolants. However, if the 
PSF is undersampled, aliasing will occur, and the PSF 
interpolation will be poor. In this situation, it is bet- 
ter to provide an oversampled PSF to Galfit (and to 
specify the amount of oversampling) , even if the data 
are undersam pled. With HST da ta this can be done us- 
ing TinyTim (jKrist fc H ook 1997j) or by combining stars. 
Galfit will take care of rebinning during the fitting. 

Note that the alternative to fitting a PSF is to fit a 
Gaussian with a small width (e.g., 0.4-0.5 pixels), which 
Galfit will convolve with the PSF. This is generally not 
advisable if a source is a pure point source because con- 
volving a narrow function with the PSF will broaden out 
the overall profile, even if slightly. The convergence can 
also be poor if the FWHM parameter starts becoming 
smaller than 0.5 pixels. However, this technique can still 
be useful to see if a source is truly resolved. 

The Background Sky The background sky is a flat 
plane with flux gradient along x and y directions. Thus 
it has a total of 3 free parameters. The pivot point for 
the sky is fixed to the geometric center (x c , y c ) of the im- 
age, calculated by (n p - lx + l)/2, where n p ; x is the number 



of pixels along one dimension. The tip and tilt are calcu- 
lated relative to that center. Because the galaxy centroid 
located at (x, y) is in general not at the geometric center 
(x c , y c ) of the image, the sky value directly beneath the 
galaxy centroid is calculated by: 



sky(£, y) = sky(x c , y c ) + (x - x c ) 



cfeky 
dx 



(y - y c ) 



dsky 
dy 
(21) 



4. THE AZIMUTHAL SHAPE FUNCTIONS 



Whereas the radial profile governs the decline of galaxy 
flux radially from a central peak, the azimuthal functions 
generate the projected shape in the x — y plane of the 
image. For instance, ellipsoidal, irregular, spiral, disky, 
and boxy shapes are all created by azimuthal functions. 
All traditional 2-D image-fitting techniques use an ellipse 
as the fundamental shape, which is obtained by stretch- 
ing the coordinate grid along one dimension compared to 
the orthogonal direction. Indeed, all azimuthal functions 
are coordinate transformations. Therefore, to change a 
shape from an ellipse into more exotic shapes, the coor- 
dinate system [r(x, y)] can be further stretched or shrunk 
radially from the peak, as a function of azimuth angle. 
This coordinate transformation preserves the functional 
form of the surface brightness profile in every direction 
because the profiles are self-similar — that is, they are 
functions of (7"/r sca i ). Thus defined, the radial profile 
parameters (e.g., r e , q, central concentration, etc.) re- 
tain their original meaning no matter the complexity of 
the azimuthal shape. 

We introduce four new ways to modify the azimuthal 
shape of a model, beginning with the traditional ellip- 
soidal model. On top of an ellipsoid, this section de- 
scribes how one can add Fourier modes, bending modes, 
and coordinate rotation functions (power law and loga- 
rithmic). Each component can be modified by any one or 
all of the azimuthal functions simultaneously, depending 
on the complexity of the galaxy one is trying to analyze. 
The next section will cover truncation functions. 

Generalized Ellipses The simplest azimuthal shape 
in Galfit is the traditional generalized ellipse. This is 
the starting point for all Galfit analysis, no matter how 
complex is the final outcome. The radial coordinate of 
the generalized ellipse is defined by: 



r(x,y)= I \x-x \ c ° +2 + 



y - yo 


Co+2\ 


^0 + 2 


q 







(22) 



Here, the ellipse axes are aligned with the coordinate 
ax es, and (xp,yo) is the cent roid of the ellipse. Defined 
by I Athanassoula et al.l (|1990D . the ellipse is called "gen- 
eral" in the sense that Co is a free parameter, which 
controls the diskiness/boxiness of the isophote. When 
Co = the isophotes are pure ellipses. With decreasing 
Co (Co < 0), the shape becomes more disky (diamond- 
like), and conversely, more boxy (rectangular) as Co in- 
creases (Co > 0) (see Figure [3]). The major axis of the 
ellipse can be oriented to any PA. Thus, there are a total 
of 4 free parameters {xo,yo, q, #pa) in the standard ellipse 
and an additional one, Co, for the generalized ellipse. 




Fourier Modes Few galaxies look like perfect el- 
lipsoids, and one can better refine the azimuthal shape 
by adding perturbations in the form of Fourier modes. 
The Fourier perturbation on a perfect ellipsoid shape is 
defined in the following way: 

r(x, y) = r (x, y) ^1 + ^ a ™ cos ( m (^ + fim))) • 

In the absence of Fourier modes in the parenthesis, the 
r o(x,y) term is the radial coordinate for a traditional 
ellipse, and 9 — arctan ((y — yo)/((x — Xo)q)) defined in 
Equation 1221 The Fourier amplitude for mode m is a m . 
Defined as such, a m is the fractional deviation in radius 
from a generalized ellipse of Equation[22l The number of 
modes N is unrestricted, and any mode can be left out. 
The "phase angle," 4> m , is the relative alignment of mode 
m relative to the PA of the generalized ellipse; the phase 
angle is 0° in the direction of the semi-major axis of the 
generalized ellipse (rather than up), increasing counter- 
clockwise. Figure [9] shows some examples of how Fourier 
modes modify a circle and an ellipse into other shapes. 

Each Fourier mode has 2 free parameters, a m and <f> m , 
and the number of modes the user can add is unre- 
stricted. However, the most useful modes are low-order 
ones (m = 1,3. ..6). We note that the m = 2 mode is 
partially degenerate with the classical axis ratio param- 
eter, for an ellipse. Therefore the use of to = 2 and q, 
together, should be largely avoided except in some situ- 
ations (e.g., peanut- looking bulges). 

The phase angles of the Fourier modes are also useful 
information to keep in mind. Modes with the following 
phase angles have the following symmetry properties: 

• Symmetry about a central point: a\ = 0, regardless 
of other mode phase and amplitude. 

• For all modes m, there is reflection symmetry at: 



<t>m = 0°j ±i^-. For to = even, this symmetry is 
about both the major and minor axes, whereas for 
to =odd, the reflection symmetry is only about the 
major axis. 



• For odd modes of to, there is additional reflection 
symmetry about the minor axis at: 



= ±90° 



An irregular galaxy has angles that are "out of phase" 
whereas regular galaxies have angles that are more "in 
phase" (i.e. reflectionally symmetric around either mi- 
nor or major axis). Therefore, it is possible to quantify 
various forms and degree of asymmetry by constructing 
indices based on the amplitude and phase angles of the 
Fourier modes. The most intuitively obvious asymmetry 
index is the to = 1 mode, which captures the lopsidcd- 
ness {Al) of a galaxy, i.e. the positioning of the bright- 
est central region relative to the fainter outer region of a 
galaxy: 

A L = \ ai \. (24) 

Asymmetric galaxies are also characterized by overall de- 
viation from an ellipse; thus, another intuitively useful 
quantity to measure is the sum of the Fourier amplitudes: 



-4, 



JV 

Ei 



(25) 



Asymmetric galaxies by definition have high Ae- How- 
ever, it is possible for galaxies with both high Ae and Al 
to be reflectionally symmetric; the degree of reflcctional 
symmetry may be an indicator for how well the galaxies 
is relaxed. Reflection asymmetry is given by the index 
A R : 

A R = Em=ovcn \ a m\ Sm 2 (tTTO^) + 

E m =odd M sin2 ( 7rTO ^) i ( 26 ) 
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where <p m is in degrees. In this formulation, the higher 
the reflectional asymmetry, the higherthe index Ar. 
Used together, these three descriptors provide highly use- 
ful ways to quantify the degree galaxies are irregular. For 
instance, high values of Ar and Al most likely imply 
high global asymmetry in the intuitive sense. Whereas 
a high value of Ae with low Ar implies high regularity, 
but large deviation from an ellipse, such as edge-on disky 
galaxies or a disky/boxy ellipticals. 

Bending Modes Bending modes allow for power- 
law-shape curvatures in the model, as opposed to spi- 
ral windings. The coordinate transformation (x, y) ==>■ 
(x',y') is obtained by only perturbing the y-axis (in a 
rotated frame) in the following way: 



v = y + a,n 



7*scalc 



(27) 



where x 1 = x, r sca i e is the scale radius of the model (i.e. 
r e g for Sersic, r s for exponential, etc.). Some examples 
of this perturbation are shown in Figure [TO] Note that 
m = 1 resembles quite closely to the axis ratio parameter, 
q. However, the m = 1 bending mode is actually a shear 
term, the effect of which is most easily seen when it op- 
erates on a purely boxy profile with Co ss 2 (Figure EJj) , 
shearing it into a more disky shape (see Figure HUtT). 
The bending modes can be modified by Fourier modes 
or diskiness/boxiness to change the higher order shape 
of the overall model. This kind of coordinate transfor- 
mation again preserves the original meaning of the radial 
profiles. Here, the object size parameter refers to the 
unstretched size, i.e. projected onto the original (x, y) 
Cartesian frame, as opposed to a length along the curva- 
ture. 

Coordinate Rotation: The Concept Sometimes 
the isophotes of a galaxy can rotate as a function of ra- 
dius, as in the case of spiral galaxies. To model spiral pat- 
terns, it is now possible to allow for coordinate rotation in 
Galfit. Coordinate rotation in Galfit means that the 
flux within circular annuli overlayed on a model rotates 
as a function of radius, i.e. 9 — f(r). The functional form 
f(r) can be fairly arbitrary but the most familiar pattern 
in nature is that of a logarithmic spiral, i.e. 9 ~ log(r). 
However, many spiral galaxies deviate from logarithmic 
winding either in the inner region, for instance due to 
the presence of a bar, or in the outer region, as might 
be due to tidal or non-relaxed features. These structures 
pose a problem when fitting galaxy images because one 
cannot simply mask out regions of non-interest when the 
goal might be to obtain the cleanest separation between 
a spiral and other embedded components. Therefore a 
pure logarithmic spiral, while useful to trace segments of 
a spiral, is often not ideal for fitting the whole galaxy, 
but ought to be modified in some ways. For this reason 
we introduce the concept of a hyperbolic-tangent (tanh) 
modification to a logarithmic or a power-law spiral. 

A pure tanh function looks like Figure 111b . showing 
that f(r) asymptotes to constant values at r — > ±oo, 
which is a highly desirable feature. As shown in Fig- 
ure [TTla. the function can be scaled, stretched, and shifted 
so that 9(r) ss at r < rj n : it is useful to model a bar- 
like feature, which, by definition, has a constant PA as 



a function of radius. A tanh function is also useful in 
the upper asymptotic limit because f(r) at r > r ou t, 
when multiplied by another function /aO"), takes on the 
form of /2(f)) an d the crosstalk within r; n is minimal, 
as shown in Figure lllfo . In short, a tanh function al- 
lows for a transition between two functions: a constant 
function at r < r- ln and another r > r out , for example a 
power-law or a logarithmic function. Moreover, the rate 
of that transition can be cleanly managed and is easy to 
interpret. For this reason a hyperbolic tangent is also a 
function of choice later on in Section [5] when we present 
the idea of a truncation function. Galfit allows for two 
types of coordinate rotation functions, the power-law spi- 
ral (a-tanh), and the logarithmic spiral (log-tanh), both 
of which are motivated empirically. We note that even 
though the logarithmic spiral is favored more in the lit- 
erature, we find that the a-tanh spiral is better able to 
capture the range of spiral behaviors found in nature be- 
cause of the one extra degree of freedom in a, which can 
simulate the behavior of the log-tanh spiral over regimes 
of interest. We therefore tend to prefer use of the a-tanh 
coordinate rotation by default. We now give an overview 
of the two types of coordinate rotation: 

Coordinate Rotation I: Power-law - Hyperbolic 
Tangent (a-tanh) The term "power law" refers 

to the fact that the pure tanh function of Figure [TTJa is 
multiplied by a function of the form ~ r a . The exact 
functional form of the rotation function is lengthy (see 
Appendix A), but the schematic functional dependence 
of the power-law spiral on the parameters is given by the 
following: 



1 / r 

2 \r out 



6{r) = 9 out tanh (r in , r out , 6* inc i, 9 S pl; rj : 

' "' (28) 

As defined, the power-law rotation starts to take hold 
beyond r — r out , and below which the tanh transition 
dominates. Figure [TT] shows a pure hyperbolic tangent 
rotation function for several different values of the pa- 
rameter rin (left)., and a combination of "bar" (n n ) pa- 
rameter and the asymptotic power- law slope a (right), 
where r is the radial coordinate system and 9 ou t is the 
rotation angle roughly at r out . The inner radius, r- lni 
is defined to be the radius where the rotation reaches 
roughly 20° . This angle corresponds fairly closely to our 
intuitive notion of bar length based on examining images, 
but is not a rigorous, physical definition. The angle 9i nc \ 
is the line-of-sight inclination of the disk, where 6*i nc i = 0° 
is face-on and 6*; nc i = 90° is perfectly edge-on. 

To motivate intuition for the free parameters used in 
the coordinate rotation definition, Figures [TH and [T3] 
show a progressive series of images for the spiral rotation 
function with different combination of parameter values. 
For instance, Figure [T2l shows a series of images of pure 
hyperbolic tanh spiral with increasing maximum rotation 
angle (f? out ), all else being held constant at the values in- 
dicated at the top. The spiral arm winding increases 
with increasing out , and the winding gets tighter, but 
the body does not expand wider because r out is fixed. 
It is also important to note that a face-on model does 
not necessarily mean that the outer-most isophotes are 
round. Rather, the ellipticity of the outer-most isophotes 
is related to the asymptotic behavior of the rotation func- 



13 



q = 1.0, a m = 0.05, <f> m = 



IT 



i i i i i i i i i i i i i i i i i i 



i i i I i i i i I i i i i i i i i i 




m=2 




































|+H+ 


+H+ 


-++++- 


HH-H- 






+H+ 


-H-H- 


-H-H- 


+H+ 






-H-H- 


+H+ 


+H+ 


1 1 1 1 



m=3 




m=4 



-- m=5 



-- m=6 






I ' i ' i I i ' i ' I ' i i i I i i i i I I I i i i i I i i i i I i i i i I i i i i iTl i i i i I i i i i I i i i i I i i i 



q = 0.5, a m = 0.5, <f> m = -45 



I i i i i I i i i i I i i i i I i i i i I 
~m=l 



i i i i I i i i i I i i i i I i i i i 
m=2 



i i i i I i i i i I i i i i I i i i i 
m=3 



































{++++- 


+H+ 


-++++- 


-++++- 






-++++- 


-++++- 


-++++- 


-++++- 






-++++- 


H+hf+HH 


-++++] 




m=4 




m=5 




m=6 




i i i i i i i i i i \i /i i i i i i i i 



Fig. 9. — Examples of Fourier modes. (Top) Low-amplitude (a m = 0.05) Fourier modes modifying a circular profile (q = 1.0) with phase 
angle <p m = 0°. (Bottom) High-amplitude (a m = 0.5) Fourier modes modifying an elliptical profile (q = 0.5) with phase angle <j> m = —45°. 



tion, which asymptotes to a constant PA beyond a radius 
of r out for a pure hyperbolic tangent (a = 0, Figure [TTT z). 
The isophotes only appear circular in the main body of 
the spiral structure when it has a large number of wind- 



ings. Figure [13] shows several other examples of barred 
and unbarred spirals, with progressively different a val- 
ues, sky inclination angle, and rotated to different sky 
position angles (@pj[)- The parameters for each grey- 
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Bending Modes 




Fig. 10. — Examples of bending modes modifying a circular profile (q = 1.0) with Co = (unless indicated otherwise). (Top row) 
Low-amplitude (a m = 0.05r™ alc ) bending modes. (Bottom row) High-amplitude (a m = 0.2r™ alc ) bending modes. Bending modes can be 
combined with Fourier modes to change the higher order shape. 



Pure Tanh Spiral (a= 0) Power Law - Tanh Spiral 




-50 50 100 150 200 -50 50 100 150 200 
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Fig. 11. — Hyperbolic tangent-power-law spiral angular rotation functions with outer spiral radius of r ou t = 100. (a) Examples of pure 
hyperbolic tangent spi rals (a = ) with different bar radii (r- ln ). (b) Examples with different bar radii and asymptotic power law a, as 
indicated. See Figures j 121 and [T3l for examples of how these parameters translate into 2-D images. 



scale figure are shown at the top and to the right of the 
corresponding (column, row). When the power-law in- 
dex a is negative, the spiral pattern can reverse course 
after reaching a maximum value (see right-most column 
of Figure 

In summary, the hyperbolic tangent power-law func- 



tion has 6 free parameters: # ut, ?~in, ?~out, a, #; nc i, and 
0^tj[. The thickness of the spiral structure is controlled 
by the axis ratio q of the ellipsoid being modified by the 
hyperbolic tangent, or by the Fourier modes that modify 
the ellipsoid. To create highly intricate and asymmetric 



15 



50 

50 



50 

50 



Pure taiih spiral: r in =o r out = so 9 il>d = o 



50 h 


50 



i i i 



T 



i i i 



I I I I I 1 1 



270 



I I I I I I I 



5-10 




l l l I l l l 



I I I 



T 



i i i 

90 




I I I I I I I 



360 




I I I I I I I 



630 




l l l I l l l 



I I I 



T 



i i i 



a Siii 



I I I I I I I 



450 



1 1 1 1 1 1 1 



720 




l l l I l l l 



50 50 



50 50 



50 50 



Fig. 12. — Examples of pure (i.e. with power law a = or without logarithmic function) hyperbolic tangent coordinate rotation modifying 
an elliptical profile with axis ratio q = 0.4. Note that all panels share the same parameters as shown up top. The spiral model has no bar. 
The numbers within each panel show the amount of total winding (units in degrees) at the spiral rotation radius of 50. Notice that outside 
r = 50, the rotation angle becomes constant, due to the rotation function being a hyperbolic tangent, thereby creating the appearance of 
a flattened disk, even though there is not a separate disk component involved in the model. 



spiral structures, Fourier modes can be used in conjunc- 
tion with coordinate rotation. 

We note that the "bar" radius (rj n ) is a mathematical 
tool. Even though the ri n term in the coordinate rota- 
tion does look like a bar when it is sufficiently positive, 
it should be regarded only as a mathematical construct 
to grant the rotation function as much flexibility as pos- 
sible. This construct can reflect reality, but it does not 
have to. For instance, mathematically, a negative r; n 
radius (Figure ITTb ) is perfectly sensible because of the 
way Equations [35] and [23 (f° r logarithmic spirals, be- 
low) are defined: a negative r m value just means that 
the spiral rotation function has a finite rotation angle 
at r = relative to the initial ellipsoid out of which it 
is constructed. When there is clearly no bar, the r- m pa- 



rameter can become quite negative; in this case, the fit is 
often indistinguishable from one where the bar radius is 
0. Furthermore, often times, one may not wish to create 
a bar and a spiral out of one smoothly continuous func- 
tion for various reasons, for instance because they may 
have different widths, the spiral may not extend into the 
center, or the spiral may start off in a ring. In these situ- 
ations, one can "detach" the bar from the spiral by using 
a truncation function (see § [5| , by instead creating a bar 
with a separate Sersic, Ferrer, or other function. When 
this is done, a "bar radius" is still useful mathematically 
in the coordinate rotation function, but it may bear no 
physical relation to the physical bar. 

Finally, we draw attention to some limitations of the 
spiral rotation formulation. While the a-tanh rotation 
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Fig. 13. — Examples of power law - hyperbolic tangent (a-tanh) coordinate rotation modifying a face-on (0j nc i = 0°) elliptical profile 
with axis ratio q = 0.4. The parameters of the rotation functions arc shown on the top and right-hand side of the diagram. The top panels 
show the spiral rotation angle as a function of radius for the panels in the same column. In the right-most column, the spiral arms reverse 
direction at r = 30 because the spiral rotation function (top-right panel) decreases in rotation angle. 



function works surprisingly well for many spiral galaxies, 
the function is smooth, so "kinks" in the spiral struc- 
ture cannot yet be modeled, even though it is possible 
to do so by allowing for "kinks" in the rotation function. 
Lastly, the spiral structure cannot wind back onto itself, 
because that would require the rotation function to be 
multi- valued. 

Coordinate Rotation II: Logarithmic - Hyper- 
bolic Tangent (log-tanh) The winding rate of 
spiral arms in late-type galaxies is often thought to be 



logarithmic with radius rather than power law. Thus, 
Galfit also allows for a logarithmic-hyperbolic tangent 
coordinate rotation function, which is defined as: 



0(r) 



sky. 



Tim ''out > C'incl > "pA > T I 



! out tanh ( 
logf— + l) /logf^ + l) 



(29) 



Like the a-tanh rotation function, the log-tanh func- 
tion has a hyperbolic tangent part that regulates the 
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Fig. 14. — Logarithmic - hyperbolic tangent spiral angular rotation functions, (a) Examples of different bar radii, where the outer 
hyperbolic spiral radius is r ou t = r- m + 10. The lower horizontal dashed line shows the rotation angle at the "bar" radius (rj n ). (b) 
Examples with different "bar" radii (r; n ) and winding-scale radii r ws , as indicated, illustrating the degree of flexibility of the spiral rotation 
rate. The rotation angle at r ou t is fixed to 150° , as shown by the upper horizontal dashed line. The left-most, black curve is close to being 
a pure logarithmic function, recasted so that at r = 0, the rotation angle 8 = 0°. 



bar length and the speed of rotation within r out . Be- 
yond r ou t the asymptotic rotation rate is that of the 
logarithm function, which has a winding scale radius 
of r ws ; the larger the winding scale radius, the tighter 
the winding. Thus, like the ct-tanh spiral, the log- 
tanh spiral rotation function also has 6 free parameters: 
6Uit, r in , r out , r ws , 9 inch and 0^1 . Note that in terms of 
capabilities, the a-tanh function can often reproduce the 
log-tanh function and more. Therefore, the a-tanh is 
probably a more useful rotation function in practice. 

Note that Galfit does not allow for a pure logarith- 
mic spiral because such a function has a negative-infinity 
rotation angle at r = 0. Therefore, in Galfit, at r = 
the rotation function reaches 9 — (Figure fl4|) . Lastly, 
it is also important to keep in mind that the meaning 
of the "bar radius," just as described in the section for 
a-tanh rotation function, is a mathematical construct. 



5. THE TRUNCATION FUNCTION 

Truncation functions allow for the possibility of cre- 
ating rings, outer profile cut-offs, dust lanes, or a com- 
posite profile in the sense that the inner region behaves 
as one function and the outer behaves as another. The 
truncation function can modify both the radial profile 
and azimuthal shape. A ring can be created by truncat- 
ing the inner region of a light profile. Likewise, when a 
galaxy has spiral arms that do not reach the center, it 
can be viewed as being truncated in the inner region. 

5.1. General Principle 

In Galfit each truncation function can modify one or 
more light profile models. Also, any number of light pro- 
files can share the same truncation function. The trunca- 
tion function in Galfit is a hyperbolic tangent function 



(see Equation [7] in Appendix B). Schematically, a trun- 
cated component is created by multiplying a radial light 
profile function, fo 7 i(x,y; ...), by one or more truncation 
functions, P m or 1 — P n (depending on whether the type 
is an inner or an outer truncation), in the following way: 

fi(x,y; ...) = f .i(x,y;x c ^,y Ctl ... qi,8pA,i) x (30) 

771 

Ar soft Qmi #PA,m) X 

n 

5 1) i Xc,7i> Vc^ni ?"break,7i5 Arsoft, m Qm ApA.n)] ■ 

The break radius, r D reak, is defined to be the location 
where the profile is 99% of the original (i.e. untruncated) 
model flux at that radius. The parameter Ar so f t is the 
softening length, so that r — rb rca k ± Ar sort is where the 
flux drops to 1% of that of an untruncated model at the 
same radius (the ± sign depends on whether the trun- 
cation is inner or outer). The inner truncation function 
(Pm) tapers a light profile in the inner regions of a light 
profile, whereas the outer truncation function (1 — P n ) 
tapers a light profile in the wings. 

The behavior of the hyperbolic tangent function is ideal 
for truncation because it asymptotes to 1 at the break 
radius r > rb rca k and at the softening radius r < r so f t , 
and vice versa for the complement function. Thus, when 
multiplied to a light profile /(r), the functional behavior 
exterior to the break radius has intuitively obvious mean- 
ings. For example, as shown in Figure 116b . if a Sersic 
function with n — 4 is truncated in the wings (shown 
in red), the core has exactly an n — 4 profile interior to 
'"break (marked with a vertical dashed line), which is a 
free parameter to fit. Likewise, an n = 4 profile trun- 
cated in the core (green) has exactly an n = 4 profile 
exterior to the outer break radius. Thus, when one sums 
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Fig. 15. — Logarithmic - hyperbolic tangent spiral (log-tanh) angular rotation examples, all face-on (0i nc i = 0°) and = 0°. The 

top-left panel shows the meaning of the rotation parameter values at the corners of each box. As with the ct-tanh spirals, the log-tanh 
spiral can be tilted and rotated to any sky projection angle, or combined with Fourier modes to produce lopsided or multi-armed spiral 
structures (not shown), and with truncation function to produce an inner ring or an outer taper. The top-left panel figure, for all practical 
purposes, is a pure logarithmic spiral with a winding scale radius r ws = 5. 



two functions of different Sersic indices n (Figure [Wb ) 
the asymptotic profiles of the wing and core retain their 
original meaning, and there is very little crosstalk out- 
side of the truncation region (denoted by vertical dashed 
lines in Figure [To]) . 

Use of the truncation functions is highly flexible. There 
can be an unrestricted number of inner and outer trunca- 
tion functions for each light profile model. Furthermore, 
multiple light profile models can share in the same trun- 
cation functions. This is useful, for instance, when trying 
to fit a dust lane (inner truncation) in a fairly edge-on 
galaxy that may affect both the bulge and the disk com- 
ponents. Just as with light profile models, the trunca- 
tion functions can be modified by Fourier modes, bending 



modes, etc., independent of the higher order modes for 
the light profile they are modifying. 

5.2. Different Variations of the Truncation Function 

Truncation models appear in many physical contexts, 
such as dust lanes, rings, spirals that do not reach the 
center, joining a spiral with a bar, or cut-off of the outer 
disk. To allow the truncation parameters to be more 
intuitive to understand given situations at hand, Gal- 
fit offers several variations. In addition to inner and 
outer truncations, truncation functions can share in the 
same parameters as the parent light profile. There are 
radial and length/height truncations, softening radius vs. 
softening length (default vs. Type 2), inclined vs. non- 
inclined (default vs. Type b) truncations, and, lastly, 
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Radius Radius 

Fig. 16. — Examples of hyperbolic truncation functions on n = 4 and n = 1 Sersic profiles, (a) A continuous n = 4 model represented 
as two truncated models of otherwise identical r e , n, and central surface brightness, with truncation radii at r = 15 and r = 20, as marked 
by the vertical dashed lines. The black curve is the sum of the inner and outer functions. This shows that, outside the truncation region, 
there is very little "crosstalk" between the inner and outer components. (6) A composite profile made up of an n = 4 nucleus truncated in 
the wings and an n = 1 truncated in the core, with truncation radii r = 10 and r = 20. Note that the hump in the summed model would 
give rise to a ring in a 2-D model. 



four different ways to normalize the flux — the most sen- 
sible choice depends on how a profile is truncated. We 
now discuss each of these variations in more detail. 

Parameter Sharing. In the most general form, each 
truncation function has its own set of free parameters: 
Xq, Ho, fbroak, Arsoft, q, and 6 PA . However, by default, the 
parameters xo,yo, q, and 9pa are tied to the light profile 
model, and are activated only when the user explicitly 
specifies a value for them. 

Radial ("radial") vs. Length ("length") or 
Height ("height") Truncations. The most use- 
ful type of truncation is one that has radial symmetry 
to first order, i.e. where it has a center, an ellipticity, 
and an axis ratio. However, in the case of a perfectly 
edge-on disk galaxy ("edgedisk" model), an additional 
type is allowed that truncates linearly in length or in 
height. For instance, a dust lane running through the 
length of the galaxy has an inner height truncation. For 
the "edgedisk" profile, Galfit also allows for a radial 
truncation, as with all other functions. The one draw- 
back to height and length truncations is that they cannot 
be modified by Fourier and higher order modes like the 
radial truncations. 

Softening Length ("radial") vs. Softening Radius 
("radial2"). Sometimes, instead of softening length 
(Ar so f t ), it is more useful for the fit parameter to be 
a softening radius (r so ft), especially when one desires to 
hold the parameter fixed. That is also allowed in Galfit 
as a Type 2 truncation function, designated, for exam- 
ple, as "radial2." The default option does not have a 
numerical suffix. 

Inclined (default, "radial") vs. Non-inclined 
("radial-b") Truncations. A spiral rotation func- 



tion is an infinitesimally thin, planar structure. Never- 
theless, it should be thought of as a 3-D structure in 
the sense that the plane of the spiral can be rotated 
through three Euler angles, not just in position angle 
on the sky. When a truncation function is modifying a 
spiral model, it is therefore sometimes useful to think 
about the truncation in the plane of the spiral model. 
When Fourier modes and radial truncations are modi- 
fying a spiral structure, the default ("radial") is for the 
modification to take place in the plane of the spiral struc- 
ture. However, there are some instances when that may 
not be ideal (e.g., a face-on spiral may actually be ellip- 
soidal). In those situations, one can choose "radial-b", 
which would allow a truncation function to modify the 
spiral structure in the plane of the sky, even though the 
spiral structure can tip and tilt as needed. 

Lastly, the truncation function can be Type 2b (i.e. 
"radial2-b") as well. 

Flux Normalization. The most intuitive flux nor- 
malization for a truncated profile is the total luminosity. 
Unfortunately, both the total luminosity and the deriva- 
tive of the free parameters with respect to the total lu- 
minosity are especially time-consuming to work out com- 
putationally and slow down the iteration process. There 
are generally no closed form analytic solutions to the 
problem. Therefore, the alternative is to allow for differ- 
ent ways to normalize a component flux. The user may 
choose whichever one is more sensible, given the situa- 
tion and the science task at hand. The default depends 
on the truncation type: 

• Inner truncation: the flux is normalized at the 
break radius. This is most sensible for a ring model 
because this radius roughly corresponds to the peak 
flux of the ring. 
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Fig. 17. — Examples of truncation functions acting on a single- component light profile of various shapes, (a) Inner truncation of a 
round profile, creating a ring, (b) The truncation function can be modified by Fourier modes, just like the light profile, (c) The truncation 
function can be offset in position relative to the light profile. (<2) The truncation function can act on a spiral model, (e) The truncation 
can tilt in the same way as the spiral. (/) The truncation function can be modified by Fourier modes while acting on a spiral model, (g) 
A round light profile is being truncated in the wing by a pentagonal (Fourier mode 5) truncation function, (h) A round light profile is 
being truncated in the inner region by a triangular function (Fourier mode 3), and in the wing by a pentagonal function, (i) A three-arm, 
lopsided, spiral light profile model is truncated in the wing by a pentagonal function, and in the inner region by a triangular function. 



• Outer truncation: flux normalized at the center. 

• Both inner and outer truncation: same as the case 
for inner truncation. 

However, there are many situations when the default 
is not desirable. Instead, the user can choose the radius 
where the flux is normalized. To be pedagogical, we ex- 
plicitly show here the normalization for just the Sersic 
function: 

• function : default (e.g., "sersic," "nuker," "king," 
etc.). Sec the details for individual functions. 



• functionl: flux normalized at the center r = 
(i.e. So). A function that is given originally by 

/orig(r) is now defined as / mod (r) 



S { orlg S . For 



'/orig(O) 

the Sersic profile (i.e. called "sersicl"), the profile 
function is redefined in the following way, written 
explicitly: 



exp 



/mod(r) = S - 



l/r 



exp [k] 



(31) 



For the Ferrer and King profiles, this normalization 
is the same as the default normalization. 
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• function2: flux parameter is the surface brightness 
at a model's native size parameter (parameter 4 of 
the light profile model). For a Sersic profile, called 
"sersic2," this means the effective radius r e . So, 

/mod(f) = S e f°" S (r\ ■ F° r example, a Sersic profile 
now has the following explicit form: 



/mod(f) = S e exp 




(32) 



For the Nuker profile this normalization is the same 
as the default normalization. 

• functionS: flux parameter is the surface brightness 
(Ebreak) at the break radius (r\, rea k)- This is the 
most useful situation when a truncation results in 
a large-scale galaxy ring, so that the surface bright- 
ness parameter corresponds closely to the peak of 
the light profile model. When the truncation is 
not concentric with the light profile model, this 
kind of normalization is not very intuitive. For 
"radial" truncation, fbreak is parameter 4, whereas 
for "radial2," rbroak is parameter 4 for outer trun- 
cation and parameter 5 for inner truncation. When 
the "sersic3" option is chosen, the rb rea k parameter 
comes automatically from the first truncation com- 
ponent with which a certain light profile model is 
associated. 

In our example of the Sersic profile, / mo d(r) — 



^brcak /or . g(rb: ^ 

has the following explicit form: 



For example, a Sersic profile now 



exp 



/modW — Sbreak" 



-K — 



1/r 



exp 



l/n 

''break \ _ 1 



(33) 



Figure [T7| demonstrates just some of the possibilities 
allowed when fitting truncations. In addition to the reg- 
ular ellipsoid shape, the higher order modes like diski- 
ness/boxiness parameters, bending modes, and Fourier 
modes can also modify the shape of the truncation func- 
tions. One can also use the truncation function on a spi- 
ral model, on models with Fourier and bending modes, 
and diskiness/boxiness models, some of which are shown 
in Figures ITm . fTTl e , 16/ , and fTTt . When a truncation 
function acts on a spiral component, it can do so cither 
in the plane of the disk ( "Type a" ) or in the plane of the 
sky ("Type b,"; e.g., "radial-b"). While the default is in 
the plane of the disk, the parameters are more intuitive 
in Type b cases when the disk is tilted and rotated. 

5.3. Caveats about using the Truncation Function 

The use of truncation functions should be carefully su- 
pervised because unexpected things can happen, such 
as the size or the concentration index of a component 
can grow without bound. This behavior is due to the 
fact that there are degeneracies between the sharpness of 
truncation and the steepness/size of the galaxy. There- 
fore, truncation functions should only be used on objects 
that clearly have truncations. 



When two functions are joined by using a truncation 
function, the crosstalk region is located in between the 
two truncation radii: it is worth bearing in mind the defi- 
nition that at the break and softening radii, the fluxes are 
99% and 1% that of the same model without truncation, 
respectively. In other words, the larger the truncation 
length, the larger the crosstalk region. Therefore, when 
one (or more) of the parameters r brC ak, ^brcak + Ar so ft, 
or r so f t is either too small (< few pixels) or larger than 
the image size, it probably indicates that profile trun- 
cation parameters are not meaningful. Rather, it more 
likely reveals the fact that there is a mismatch between 
the light profile model and the actual galaxy profile. 

6. INTERPRETATION, PARAMETER DEGENERACIES, 
UNIQUENESS, LOCAL MINIMA, AND ERROR ANALYSIS 

Now that we have introduced several ways to modify an 
ellipse into more exotic shapes, a natural question to ask 
is how unique or robust are the modifications. A single- 
component ellipsoid fit can often be used to quantify the 
global average profile of galaxies. However, beyond that, 
decisions about what procedure to use get to be more 
complicated. On the one hand, the science goal might 
call for fitting detailed structures inside a galaxy (e.g., a 
bulge, bar, nuclear star cluster, etc.). On the other hand, 
doing so raises concerns about parameter degeneracies, 
uniqueness, and local minima solutions when the analysis 
becomes complex. It is therefore useful to consider in 
some depth what causes degeneracies and the different 
contexts in which they appear. Doing so allows for better 
understanding for how to deal with them and how to 
properly interpret results from complex analysis. For, 
not all complex analyses are more suspect, nor are all 
simple analyses more robust. 

The term "degeneracy" has a specific mathematical 
connotation, namely the relation of a + b = c is degen- 
erate in a and b for a constant value of c. In the galaxy 
fitting literature, "degeneracy" is often more loosely used 
to also refer to "non-unique" or "local minimum" so- 
lutions (e.g., a fit oriented at 90° from the best orien- 
tation), or strong "parameter correlation" (e.g., sky is 
anti-correlated with the Sersic index n). We will mostly 
not make the subtle distinctions here and instead will use 
the term "parameter degeneracy" generically to refer to 
all such situations. 

However, when fitting galaxies, it is more important to 
distinguish between the aforementioned real degeneracies 
from "pseudo" ones. Real degeneracies refer to correlated 
parameters, local minima, and mathematically degener- 
ate solutions. By contrast, "pseudo" degeneracies involve 
convergence issues when an algorithm is used beyond its 
technical limits, or when users provide bad input model 
priors to fit the data. They may have nothing to do 
with real degeneracies, yet the behavior of convergence 
may seem to suggest otherwise. Whereas problems with 
real degeneracies are often resolvable by using full spatial 
information of 2-D images, pseudo-degeneracy problems 
are solved through experience and by using sound scien- 
tific or technical judgment, as we elaborate further. 

In this section, we discuss how most of the parame- 
ter degeneracy problems are avoidable with proper input 
priors and proper fitting supervision, even when large 
numbers of free parameters are involved. We also dis- 
cuss why, contrary to popular notions, when it comes to 
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avoiding model degeneracy and local minima, it is not 
sufficient to only choose a model that is the simplest. 
Rather, it is a judicious combination of simplicity and 
realism that make for the most robust solutions. Lastly, 
these discussions are intimately connected to the issue 
of error analysis because error measurements are nearly 
always dominated by systematic issues rather than pho- 
ton noise in galaxy fitting. We therefore discuss why it is 
more important to quantify mo del- dependent systematic 
errors rather than to rely on statistical estimates. 

We note that the discussions below are mostly based 
on experience, which we present using practical exam- 
ples rather than to show using rigorous proof. Carrying 
out a rigorous proof is not only beyond the scope of this 
study, but it is nearly impossible to do in a general man- 
ner because different scientific applications have different 
sensitivities to different types of degeneracies. We are 
also aware that presenting a full discussion of degeneracy 
issues lends credence to the common notion that multi- 
component analysis is dangerously complex. However, 
the reality is not nearly so dire when one has a proper 
understanding of the underlying issues and causes. 

6.1. True Numerical Degeneracies Caused by Correlated 
or Non-Unique Parameters 

There are well-known situations when different param- 
eters in one or more functions are capable of modeling 
the same profile behavior. This scenario is the one most 
commonly referred to in generic discussions about model 
degeneracies. For instance, very large Sersic index val- 
ues (n > 4) have highly extended wings, the presence of 
which is non-unique with the sky parameter. A high n, 
caused by profile mismatch or poor model prior, can often 
suppress the sky estimate. It is therefore advisable to es- 
timate the sky independent of the fit, and to hold it fixed 
to the best estimate. As a second example, in the Nuker 
profile (Equation 1 17p. there are three parameters (a, /3, 
and 7) that control the inner/outer slopes and sharpness 
of the bending (Figure [7]) . When the break radius rt of 
a Nuker profile is sufficiently small and profile mismatch 
sufficiently large, model discrimination relies entirely on 
the power law . Because there are numerous ways to 

yield a specific value for in the model, it leads to a 
degenerate situation involving three parameters. 

As another example, a low-amplitude second Fourier 
mode and the first bending mode (shear) can both be 
degenerate with the axis ratio q parameter of an ellipse, 
therefore they should not be used together except in 
obvious situations where doing so is useful. Lastly, in 
the spiral rotation function, the periodicity of the rota- 
tion function can sometimes be a source of "degeneracy." 
Multiple windings can approximate a smooth continuous 
model, whether or not there is a spiral structure present. 
For instance, a classical Sersic ellipsoid can be simulated 
by a spiral model with a very large out . While the fit 
is not good and easy to diagnose by an end user, it is 
nevertheless a numerically allowed solution. 

The above situations are not meant to be a complete 
laundry list, but they are the most common situations. 
In complex analysis, one always needs to be circumspect 
about the potential hazards of mixing and matching dif- 
ferent functions whose parameters produce similar profile 
behaviors. Even though Galfit allows for a great deal of 



flexibility in the analysis, it is ultimately up to the user 
to decide on what to allow, based on the goals of the 
science, and to understand when potentially degenerate 
parameters may be used effectively. 

The above discussion may also seem to imply degen- 
eracies or non-uniqueness are too numerous for complex 
analysis to be practical or reliable. That notion is only 
true when it is not possible to verify the results of a fit 
and to try out other solutions. Such a scenario is more 
common for large scale galaxy surveys, in which auto- 
mated, detailed, analysis is admittedly quite difficult to 
conduct sensibly. However, even in those scenarios, there 
are many situations where mutually coupled parameters 
do not affect the other main parameters of scientific inter- 
est: degeneracies in the Fourier modes often do not have 
any bearing on the total luminosity or size of a compo- 
nent. Moreover, when an analysis is done manually, it 
is reassuring that the problems are almost always easy 
to recognize and remedy when they do happen, even by 
simple inspection of the model and residual images. 

6.2. Pseudo-Degeneracies Caused by Technical 
Conditions (e.g., Model Profile Resolution, 
Parameter Boundaries) 

Occasionally, what appears to be numerical degener- 
acy problems may be caused by someone using a code 
outside the algorithm's physical capabilities. As such, 
it is a pseudo-degeneracy. Different algorithms have dif- 
ferent limitations that affect convergence, whether the 
code is gradient descent, Metropolis, or otherwise. This 
situation may appear like parameter degeneracy because 
restarting the fit does indeed yield a different solution, 
but in fact the code may be hamstrung in its conver- 
gence. For example, gradient descent algorithms require 
the calculation of a gradient, and thus can run into prob- 
lems when the gradient cannot be calc ulated properly. I n 
simulated annealing algorithms (e.g.. lPress et al.lll992f> . 
parameter boundaries and annealing speed control the al- 
gorithmic behavior: anneal too quickly, the solution may 
settle into a local minimum. To search larger parameter 
spaces requires longer annealing times. 

While all algorithms have conditions under which they 
perform poorly, pseudo-degeneracies can always be rec- 
ognized and mitigated. Galfit is based on a Levenberg- 
Marquardt subroutine that performs the least-squares 
minimization. In part a gradient descent algorithm, the 
convergence behavior is affected by the calculation of gra- 
dient images that determines the direction of steepest \ 2 
descent. When the gradient images are problematic, they 
affect the convergence to a proper solution. There are 
three main problematic situations. The first, and most 
common, occurs when a model becomes extremely com- 
pact (FWHM < 0.5 pixel), so that the profile gradient 
cannot be resolved: all the gradient information in the 
model fits into a single pixel. This situation mostly arises 
when working with high-contrast imaging data, such as 
quasar host galaxy decomposition, when one of the sub- 
components may be used to reduce the strong residuals 
caused by a PSF mismatch. A similar situation arises 
when a model is very thin (axis ratio q < 0.05) and the 
object is compact; here, the gradient does not exist along 
one spatial direction because of a lack of pixel resolution. 
Another rare example occurs when the inclination angle 
of a spiral rotation component is close to perfectly face-on 



23 



(#inci —> in Equations [28l and |29|) . when the derivative 
image for the inclination parameter approaches zero. 

Another abnormal numerical behavior may occur when 
one places parameter constraints on a model to prevent 
some parameters from wandering too far from their ini- 
tial values. Doing so may cause poor convergence by 
forcing the solution into a tight "corner," when the best 
solution is somewhere beyond it. A typical situation is 
where there are other sources in the image that are not 
masked or fitted by models, but that are sufficiently lu- 
minous to influence the fit of the target of interest. In 
this situation, no amount of effort will produce a sensi- 
ble solution, because the best solution is outside of the 
parameter boundaries, even though the desired solution 
may be within. Pseudo-degeneracies occur in this situ- 
ation both because there is an abnormal condition im- 
posed and because the input prior for the model is poor. 

While technical issues with code operation add a layer 
of complexity to image analysis, in practice the majority 
of situations one encounters are straightforward to recog- 
nize by observing when the parameters take on extremely 
large or small values. However, clearly recognizing the 
problem as being pseudo degeneracies is the key. Once 
diagnosed, these situations are easy to guard against, 
by holding those parameters fixed when they go below 
certain values. In practice, technical issues are not prob- 
lematic even when Ga lfit is used for automated analysis 
(jHaussler et al.ll2007t) 7 . 

6.3. Pseudo-Degeneracies Caused by Bad Input Model 

Priors 

One of the most common causes of degeneracy prob- 
lems in galaxy fitting analysis comes from using input 
priors that are not well suited to the data. The most 
common "input priors" involve the choice of the type or 
the number of components in a model 8 . Input priors are 
ideal when the number of components of a model used 
in a fit matches the number of luminous components in 
a galaxy. However, often times one may choose to use ei- 
ther fewer or more components than needed by the data. 

A common example where the input prior is bad is 
when one uses fewer model components than called for 
by the data. Two of the main reasons for doing so are to 
reduce the number of components/free parameters, or to 
allow automated analysis, where it is not yet possible to 
tailor fits to individual galaxies. This approach is often 
an intentional course of action taken by many studies, 
especially when it comes to automating two component 
analysis; the goal, ostensibly, is to decompose a galaxy 
into bulge and disk components. Seemingly reasonable 
and justifiable on the notion of reducing the potential for 
degeneracy, the approach is generally regarded by most 
people to be a positive attribute, rather than a source 
of problem itself. Yet, that intuitive notion conflicts 
with the basic principle of how least squares algorithms 
work, and leads to perhaps the most common causes of 
(pseudo-) degeneracy problems cautioned by literature. 

7 While these conditions can always be anticipated in advance, 
implementing a solution in the code is more tricky, because the act 
of doing so may also induce other convergence difficulties. This 
leads to a false sense of security about the robustness of a solution. 

8 An input prior does not refer to the accuracy of the initial 
parameter guesses. 



To understand why using fewer components than nec- 
essary is bad, it is important to appreciate that galaxy fit- 
ting analysis is fundamentally flux weighted. Thus, when 
a luminous structure is not accounted for, other subcom- 
ponents try to compensate, however imperfectly, for its 
presence. For instance, one may use a two-component 
model fit to a galaxy that has a bulge, disk, and bar. 
Doing so may have several different outcomes. One solu- 
tion is where one component is a sum of (disk+bar) while 
the other is the bulge. Another can be (bulge+bar) and 
disk, or perhaps a compromise (e.g., bulge -I- 0.7 bar; disk 
+ 0.3 bar). Which scenario occurs depends on the rela- 
tive contrast (i.e. flux weighting) of the bar to the bulge 
and disk, and potentially on the initial parameters of the 
three components; small perturbations may "bump" the 
solution out from one minimum into another. It is quite 
possible for there to be a "global minimum" solution to 
this problem. However, when the most meaningful solu- 
tion, physically, is simply dis-allowed by the input prior, 
a globally minimum x 2 cannot lend much credence to the 
reality of the model components. 

An input model prior might also be bad if the model 
involves using more subcomponents than inherently 
present in a galaxy. In this situation, the results de- 
pend strongly on the degree of profile mismatch between 
the model function and the data. If there is significant 
mismatch, all the components cooperate to reduce the 
residuals. For instance, it is always possible to fit mul- 
tiple exponential models to a single-component de Vau- 
couleurs profile. If the goal is to obtain the total flux, 
the sum would do a better job than using a single expo- 
nential. However, individually, the structural parameters 
may not have much physical meaning. 

Another example involving model prior is in the area 
of high-contrast imaging, where the goal is to deblend 
a central, unresolved, point source from a diffuse under- 
lying object (e.g., quasar and host galaxy). To do so 
reliably requires an accurate PSF model for the unre- 
solved source, or else the residuals may overwhelm the 
extended object, causing unreliable fits. Here, the prior 
is the PSF model. Quantifying how the prior affects the 
fitting results involves trying out different PSFs, or to in- 
clude extra components to account for the PSF residuals, 
depending on the science goal. 

These examples illustrate some of the most common 
situations where the reliability of a fit depends less on the 
number of free parameters, and more on having a proper 
model to describe the data. Beyond a single-component 
analysis, the need to make such a decision means that 
it will be difficult to automate highly detailed decompo- 
sitions of galaxies. However, while multi-subcomponent 
fitting is difficult to automate, it is reassuring that mak- 
ing a wise decision, interactively, is often not particularly 
difficult when a science goal is clearly defined. More- 
over, for galaxy surveys, where the goal is to fit single- 
component profiles to galaxies, multi- object decomposi- 
tion is quite feasible to automate (e.g., M . Barden et al. 
2009, in preparation; lHaussler et aLl l2007). 

In summary, pseudo-degeneracy conditions exist be- 
cause least-squares fitting fundamentally involves flux 
weighting: when luminous flux distributions are present 
in an image, the models are attracted toward them to 
reduce the residuals. Therefore, when all components 
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are not properly modeled, the result may be tricky to 
interpret not because of potential for model degenera- 
cies, but that the solution may have no physical meaning 
even if there is a global minimum. The solution is to in- 
crease the complexity of the analysis progressively until 
all luminous components are properly accounted. This 
process does not imply, however, that it is necessary to 
account for every component inside a galaxy for the so- 
lution to have any meaning, only that components of 
similar flux ratios ought to be simultaneously accounted 
in detailed analysis; with a few exceptions (e.g. locally 
dominant features like nuclear star cluster, nuclear ring), 
components with low fluxes generally do not significantly 
disturb the parameters of the much more luminous sub- 
components. 

6.4. Parameter Degeneracies Can be Broken by Spatial 
Information in 2-D 

One of the most common notions regarding fitting de- 
generacy is that the more free parameters there are the 
greater is the potential for degeneracy problems. How- 
ever, the sheer number of parameters is often not itself an 
indication of a potential for parameter crosstalk. Con- 
sider, for instance, that it is equally robust to fit thou- 
sands of well-separated stars as it is to fit an isolated 
one. The same is true for galaxies, even though they are 
considerably more exten ded and may ov e rlap: in large- 
scale image simulations, lHaussler et al.l (|2007l ) studied 
automated batch analysis of galaxies using one Sersic 
profile per galaxy. They find that simultaneously fitting 
overlapping or neighboring objects using multiple com- 
ponents (often 3-10 Sersic models at a time) recovers the 
input simulated parameters more accurately than fitting 
a galaxy singly while masking out the neighbors. 

Indeed, spatially well-localized sources, like a bar, ring, 
or off-nuclear star clusters, are virtually free from degen- 
eracies caused by crosstalk with other components. A 
galaxy bar is well determined because it is more elon- 
gated, has a flatter radial profile, and is more sharply de- 
fined than the surrounding bulge and disk components, 
despite being embedded within. Compact objects that 
are off-centered may also be well determined if the rest 
of the galaxy can be modeled accurately. Contrary to 
notions that more model components lead to greater de- 
generacy, it is important to consider the qualitative as- 
pects of those components: not accounting for strong 
features explicitly can yield a less reliable and less physi- 
cally meaningful fit because the solution is a compromise 
between the different subcomponents. 

6.5. Measurement Uncertainties, Parameter 
Correlation, and Parameter Degeneracies 

The issue of parameter degeneracies closely ties into 
the topic of measurement uncertainties, especially when 
the result of the analysis may depend on the input model 
in fitting galaxies. When the model fits the data perfectly 
(i.e. the residuals are only due to Poisson noise) it is 
possible to infer parameter uncertainties from the covari- 
ance matrix of free parameters, which is produced during 
least-squares minimization by the Levenberg-Marquardt 
algorithm. In galaxy fitting, ideal situations are often 
not realized because the differences between the data and 
the model profile involve not only random (e.g., Poisson) 
sources, but also systematics from non-stochastic (e.g., 



profile function or shape mismatch, neighboring galax- 
ies, etc.), and stochastic factors (overall smoothness, for 
instance due to star clusters) . The one exception is under 
low signal-to- noise (S/N) situations, when Poisson noise 
exceeds model imperfection. In most other situations, 
non-random factors dominate the residuals, causing un- 
certainties inferred from covariance matrices to be un- 
derestimated. Therefore, it is frequently not very mean- 
ingful in galaxy fitting to cite measurement uncertainties 
for the fitting parameters in the traditional sense. 

One way to quantify uncertainties, possible in large 
galaxy surveys, is to allow the scatter of the data points in 
physical relations (e.g., radius vs. luminosity, luminosity 
vs. mctallicity, etc.) to articulate the overall uncertainty 
of the measurements, even if individual errors could not 
be easily obtained. Such a scatter inherently involves a 
convolution of several error sources: the intrinsic scatter 
present in a physical relation, Poisson measurement er- 
ror, stochastic and non-stochastic systematic errors due 
to model imperfections. Intrinsic scatter, being a fact 
of nature, remains present in physical relationships even 
should the data have infinite S/N, and even if the models 
are perfect fits to the data. Intrinsic scatter is often a sci- 
entifically interesting quantity, but it is difficult to differ- 
entiate from scatter caused by systematic and stochastic 
errors, which do not vanish given infinite S/N. 

In the absence of large galaxy surveys, it is then impor- 
tant to quantify stochastic and non-stochastic systematic 
errors for individual objects. Some example situations 
incl ude the black hole mass vs. galaxy relation stud- 
ies (iKormendv fc Richstonej[l995t iGebhardt et al.l 120001 
20001) and the fundamental plane 
1987) 
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In general it is very difficult to pin-point all the causes 
of non-stochastic systematic errors in an analysis, and to 
quantify their magnitude. However, one common cause 
is profile model mismatch: to the extent that one does 
not know the intrinsic model of a galaxy a priori, the un- 
certainty in measuring the parameters is wedded to one's 
assumptions about the model. Therefore, the process of 
quantifying systematic, mo del- dependent errors involves 
exploring the degree of parameter coupling, by trying 
out different models and seeing how the parameters of 
key scientific interest change. Another source of system- 
atic error is due to comparing results from different al- 
gorithms. In this scenario, the most sensible practice is 
therefore to only compare parameters that are derived 
using the same fitting technique (rather than 1-D vs. 2- 
D), and using the same pixel and flux weighting scheme 
(instead of Poisson vs. non-Poisson) during analysis. 

In contrast, stochastic errors arising from general non- 
smoothness of a galaxy profile are caused by, for example, 
star forming patches, dust lanes, etc.. Existing on small 
scales and widely dispersed, non-smoothness cannot be 
easily identified and modeled in a practical manner using 
multiple components. Even if it is possible to do so, 
whether they ought to be fitted explicitly, masked, or not 
treated at all, falls under the purview of the science goal. 
Stochastic fluctuations often influence the analysis in a 
manner analogous to having large correlated noise in the 
data. If the fluctuations can be quantified, one possible 
solution is to include them in the fit as a variance term 
of x 2 (Equation [T]). To estimate the fluctuations requires 
first obtaining a smooth underlying model, which is not 
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always easy to do if galaxies have steep and / or irregular 
profiles. 

While it is generally difficult to disentangle stochastic 
from non-stochastic sources of systematic errors, there 
also do not seem to be obvious benefits for doing so from 
a scientific standpoint. For most applications, one should 
only be interested in the overall magnitude of the sys- 
tematic errors in a collective sense. One way forward 
is therefore to understand which parameters are most 
strongly coupled, then compare the results of different 
solutions judging by which ones are physically plausible. 
For instance, one common interest in bulge-to-disk de- 
composition is to quantify the uncertainty of the Sersic 
index n. We know that the Sersic index n takes on a 
large value when a profile has both a steep core and an 
extended wing (see Figure[3]). Therefore, quantifying sys- 
tematic errors in measuring the Sersic n might involve 
masking or fitting nuclear sources/neighboring contami- 
nation, trying out different PSFs, or fitting the disk by 
allowing for different disk Sersic index values. Properly 
judging the causes of systematic errors and accounting 
for them often would lead to more natural fits and more 
sensible parameter values, without the need to hold cer- 
tain parameters fixed to preconceived answers. 

In exploring the parameter space as described, there 
is often a concern that parameter degeneracies are too 
numerous or problematic to understand, which brings 
the discussion back full circle. As discussed in previ- 
ous sections, when the cause of parameter degeneracy is 
properly identified, and when the model priors are well 
conceived, our experience has been that spatial informa- 
tion in 2-D can often effectively break many potential 
degeneracies between the components. Even when the 
size, luminosity, and central concentration, of the dif- 
ferent components correlate they often interact in fairly 
superficial ways, and do not dramatically change what 
the model components represent physically. However, 
in situations where cross-talk is significant and there is 
no reason to prefer one solution over another (when the 
input prior is befitting), then differences in the answer 
speak to the degree of the parameter uncertainty that is 
of key interest to quantify, rather than to avoid, because 
ultimately the models are empirically motivated. 

In summary, to the extent that the results may de- 
pend on model assumptions, parameter exploration is 
the only viable way to quantify true measurement errors 
in the fit parameters. Thus, when used properly, param- 
eter coupling/degeneracy, rather than complicating the 
interpretation, offers a deeper insight into the reliabil- 
ity of the overall analysis. We illustrate the above ideas 
more explicitly in the following examples. 

7. EXAMPLES OF DETAILED DECOMPOSITION 

To demonstrate how to use the new features to ex- 
tract complex structures, we analyze five galaxies that 
are well resolved: IC 4710, an edge-on disk galaxy, Arp 
147, M51, and NGC 289. These galaxies are chosen be- 
cause they represent examples where traditional analysis 
using perfectly ellipsoid models tend to leave some ques- 
tion as to what is physically being measured and to the 
robustness of the photometry and decomposition. The 
primary purpose here is to illustrate the basic building 
blocks of galaxy morphology, not to address what are the 



most "scientifically interesting" or useful applications — 
the scope of which is far too broad to address. As such, 
each individual example is not intended to necessarily be 
"interesting" in its own right. For instance, while param- 
eterizing a ring galaxy like Arp 147 may not itself be too 
worthwhile scientifically, the concept has other relevance 
to deblending Einstein rings from lensing galaxies in the 
image plane of strong gravitational lenses, or separat- 
ing a ring from a bulge, disk, and bar in spiral galaxies. 
Indeed, these are heuristic examples meant to generate 
new ideas for potentially interesting applications, and to 
illustrate the dynamic range of capabilities in our new 
approach. 

Another goal of this section is to illustrate two seem- 
ingly contradictory notions when it comes to galaxy mor- 
phology analysis: 

• Sometimes it is not necessary to perform "full- 
blown" analysis, including spiral structures, 
Fourier modes, rings, etc.. The detailed analysis 
below will show when it is not necessary to utilize 
the full machinery in order to meet the science 
requirements, such as when the interest is to only 
quantify global properties. However... 

• Sometimes it is necessary to perform full-blown 
analysis. In situations where detailed decomposi- 
tion matters (e.g., quantifying bulge-disk-bar frac- 
tions) the most reliable analysis is to make full use 
of the machinery available. 

Indeed, the availability of new tools does not in any 
way invalidate or weaken the conclusions of hundreds of 
studies that came before this one — quite the contrary. 
Rather, the main message is that given the new capabil- 
ities, it is more important now than ever to weigh the 
relative benefits of sophistication against the drawback 
of increased difficulty and time, whereas no such options 
existed before. 

7.1. IC 4710 

IC 4710 is an SB(s)m galaxy, which has a bar- like 
feature in the middle of a roundish outer structure, as 
shown in the i?-band image of Figure II 81 which comes 
from the CINGS (Carnegie-Irvine Nearby Galaxy Sur- 
vey) project 9 . Prior to the analysis, we masked out the 
stars using the SExtractor software (|Bertin fc Arnoutsl 
1996). We analyze this galaxy using, for comparison, 
both one- and two-component regular and higher order 
models with Fourier modes, shown in Figures [FSb -i The 
best-fit parameters are given in Table [TJ which illustrates 
three different sets of analysis parameters: best fit using 
two components (Figure [18b ). a model using just the tra- 
ditional ellipsoid component (Figure [LSI?), and the same 
single-component model with Fourier modes added (Fig- 
ure \T8b ) . Figure IT8t shows the radial surface brightness 
profile of the data and the individual subcomponents of 
the best model. 

There are several points to understand from comparing 
detailed and simple analyses. The best-fitting ellipsoid 
model (Figure \TEb ) is oriented more parallel to the bar- 
like, higher surface brightness component than the lower 

9 http : //users . obs . carnegiescience . edu/lho/projects/CINGS/CINGS .html 
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Fig. 18. — Detailed analysis of IC 4710. (a) Original data. (6) Best single-component Sersic profile fit with Fourier modes m = 1 
to m = 10. (c) Best two-component Sersic profile fit each with Fourier modes, corresponding to the parameters shown in Table [T] (d) 
Best- fit residuals, (e) Component 1 of 2 in the best-fit model of Panel (c). (J) Component 2 of 2 in the best-fit model, (g) A traditional 
single-component ellipsoid fit. (h) Residuals from the model in Panel (g). (i) 1-D surface brightness profile. The individual components are 
shown as dashed lines, and the solid line coursing through the data is the sum of the two components. The lower panel shows the residuals 
of data — model. 
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Table 1. IC 4710 Fitting Results 





# 


— sersic — 


Ax ["] 


As/ ["] 




mag 




r. ["] 


n 




<? 


0PA 


deg 




Comments 






fourier 


— 


mode: ampl. &; ph 


ise [deg] 


mode: ampl. & phase [deg] 


mode: ampl. 


k phase 


deg 






Best fit 


1 


— sersic — 


0.00 




0.00 




13.71 




48.92 


0.55 




0.32 


-63.36 










0.16 




0.08 




0.00 




0.17 


0.00 




0.00 




0.07 




Inner 




fourier 




1 


0.16 


1 


-97.40 


3 


0.17 


3: -17.95 


1 


0.06 


4: 17.67 




component 








1 


0.00 




1: 1.18 


3 


0.00 


3: 0.21 


1 


0.00 


4: 


0.32 








fourier 




5 


0.05 




5: 18.37 


6: 


-0.06 


6: 13.22 


7 


0.03 


7: - 


1.20 












5 


0.00 




5: 0.34 


6 


0.00 


6: 0.23 


7 


0.00 


7: 


0.63 








fourier 




8 


0.05 




8: 10.92 


9 


0.01 


9: 4.55 


10 


0.03 


10: - 


6.68 














8 


0.00 




8: 0.17 


9 


0.00 


9: 1.63 


10 


0.00 


10: 


0.29 




Outer 


2 


— sersic — 


1.97 




26.28 




12.49 




57.24 


0.37 




0.90 


41.38 




Component 






0.24 




0.19 




0.00 




0.09 


0.00 




0.00 




0.85 








fourier 




1: 


-0.31 


1 


-39.25 


3 


0.03 


3: 55.46 


4 


0.03 


4: -27.65 












1 


0.00 




1: 0.86 


3 


0.00 


3: 1.31 


1 


0.00 


1: 


0.89 








fourier 




5 


0.04 


5 


-15.73 


6 


0.02 


6: -12.47 


7 


0.01 


7: 16.56 












5 


0.00 




5: 0.87 


6 


0.00 


6: 1.15 


7 


0.00 


7: 


1.74 








fourier 




8: 


-0.03 


8 


-13.49 


9 


0.01 


9: -19.63 


10 


0.02 


10: -15.64 












8 


0.00 




8: 0.86 


9 


0.00 


9: 0.91 


10 


0.00 


10: 


0.92 








merit 


x 2 = 


167438.77 


N dot = 127966 




JVfrco = 53 




Xl^ 


1.31 








Single 


1 


— sersic — 


0.00 




0.00 




12.15 




60.37 


0.69 




0.82 


-63.51 




component 






0.06 




0.05 




0.00 




0.12 


0.00 




0.00 




0.31 








merit 


x 2 = 


247304.81 


AW = 128009 




iVftee = 10 




xi = 


1.93 








Single 


1 


— sersic — 


0.00 




0.00 




12.15 




59.00 


0.69 




0.83 


-( 


4.43 




component 






0.12 




0.09 




0.00 




0.09 


0.00 




0.00 




0.24 




with 




fourier 




1: 


-0.05 




1: 72.02 


3: 


-0.07 


3: 27.55 


4: 


-0.02 


4: - 


7.37 




Fourier 








1 


0.00 




1: 2.39 


3 


0.00 


3: 0.31 


4 


0.00 


4: 


0.48 




modes 




fourier 




5 


0.02 




5: 5.54 


6: 


-0.01 


6: -7.72 


7: 


-0.02 


7: 15.40 












5 


0.00 




5: 0.54 


6 


0.00 


6: 0.80 


7 


0.00 


7: 


0.33 








fourier 




8 


0.01 




8: 0.28 


9 


0.01 


9: 3.16 


10: 


-0.02 


10: 


1.12 












8 


0.00 




8: 0.50 


9 


0.00 


9: 0.63 


10 


0.00 


10: 


0.25 








merit 


x 2 = 


235140.44 


Af do f = 127991 




JVfrcc = 28 




Xt- 


1.84 









Note. — Best-fitting parameters for IC 4710. The meaning of the object parameters is shown at the top for each model component. The 
statistical uncertainties for each model component, based on the covariance matrix of the fit, arc shown in the row underneath the best-fitting 
model parameters. Systematic uncertainties due to imperfect model-data match are typically 1%— 10% for the fluxes, 10%-20% for the sizes, and 
20%— 30% for the Sersic index. For the Fourier modes, the phase angle is relative to the major axis of the light profile component. Note that the sky 
parameters are not shown. The "Best /if parameters (top section) correspond to Panel (c) in Figurc [l8l "Single component" parameters (middle 
section) correspond to Panel (g), and u Single component with Fourier modes" parameters (bottom section) correspond to Panel (b). 
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surface brightness body; however, the ellipsoid model is 
much broader than the bar (Figure IT8l s) . This happens 
because a single-component fit is a compromise between 
the various subcomponents of a galaxy, and, as such, it 
reflects neither one perfectly. Allowing the azimuthal 
shape to change by adding 9 Fourier modes results in 
a shape shown in Figure \Wb . Note that because the 
profile is restricted to having a Sersic functional form in 
every direction radially from the peak, the shape does 
not have complete freedom to take on any shape, as op- 
posed to a shapelet or wavelet-type Fourier inversion: it 
is merely a higher order perturbation of the best-fitting 
ellipse. Indeed, in comparing single-component fit pa- 
rameters in Table [T]for the two models, the main Sersic 
structural parameters hardly budged, despite the Fourier 
model having 18 more free parameters. Therefore, the 
marginal returns in using more free parameters is negli- 
gible in this situation when it comes to the main Sersic 
structural parameters. However, if the scientific interest 
is to quantify the global symmetry, then the higher order 
modes are of interest. 

Another point of interest is how higher order models 
affect the accuracy of the global photometry It is nat- 
ural to expect when a model is unrealistic for a galaxy 
that the photometry is also unreliable. In Figure [T8l c, 
it is evident that a two-component model is more ap- 
propriate than the single-component fits of Figure 118b 
and 118b . However, when the flux of the two-component 
model is summed, one finds that the difference with the 
single-component fits is only 0.03 mag. This and subse- 
quent examples illustrate empirically that the process of 
least-squares minimization using even naive ellipsoids is 
often capable of providing accurate photometry to within 
0.1 to 0.2 mag, even if the galaxy shape departs from el- 
lipsoid models quite drastically. 

Lastly, Figures [TSle andfTST demonstrate that it is quite 
feasible to unambiguously disentangle embedded com- 
ponents that have different shapes, using higher order 
Fourier modes. Despite there being a large number of 
parameters, it is visually clear based on Figures [T51s and 
I18lf that parameter degeneracy is not an issue, because 
the shapes of the components are quite different. In part, 
this is possible because of how Fourier modes are imple- 
mented in Galfit: the profile function is preserved in 
every direction radially from the peak, even in situations 
where the shape is irregular, as in Figure [TSl s. 

7.2. GEMS Edge-on Galaxy 

This edge-on galaxy (Figure Q1J1 Table [3]) comes from 
the G EMS (Galaxy Evolution from Morphology and 
SED, (|Rix et al J 120041 )) project, which is an HST imag- 
ing survey of the Chandra Deep Field-South. Belying 
a benign morphological appearance is a dust lane (Fig- 
ure fT9l s) that courses through the center, complicating 
the traditional ellipsoid fitting technique. 

The analysis of even this simple object can be quite 
involved. The best-fitting model involves three compo- 
nents: a fairly compact bulge, an edge-on disk compo- 
nent, and an puffy stellar halo enveloping both. Since the 
halo component is more luminous than the bulge compo- 
nent, a two-component model fit would naturally ascribe 
the halo component to the bulge, despite there being a 
distinctly rounder component at the center. Like the 
previous example, each of the three components (Fig- 



ure I19lf -/i) is modified by Fourier modes. Furthermore, 
the best fit includes an actual model for the dust lane 
(component 4, Table [2]). The dust lane is modeled by an 
inner truncation function as discussed in Section O 

A truncation model is shown as a model "component" 
in the fit; it is unique because it is not a light profile 
model, and one cannot generate an image to see what it 
looks like. Instead, its influence is to be seen on all the 
light profile models (i.e. components 1-3; Figure \TW-h), 
where it reduces the light by the same fraction for all 
components. In every other way, the truncation function 
behaves exactly like a light profile model: it can have its 
own centroid (or not), and it can be modified by Fourier 
modes, as shown in Tabled The benefit of using a single 
truncation model for all three light profile models is not 
only to reduce the degrees of freedom, but it is also phys- 
ically motivated because foreground dust attenuates all 
background light sources by an equal fractional amount. 
Nevertheless, if desired, it is also possible to allow each 
component to be attenuated differently. 

This example also demonstrates how the result of the 
analysis depends on the input prior of the model. In the 
fit using traditional ellipsoid parameters, a mask is used 
to minimize the effect of the dust on the analysis. Yet, 
the effects cannot be completely removed. As shown in 
Table [2J the inclusion of the truncation model can sig- 
nificantly affect the structural parameters: the surface 
brightnesses can differ by 0.8 mag arcsec -2 , and the sizes 
by 10%-20%, even in this seemingly uncomplicated situa- 
tion. Moreover, the differences far surpass the statistical 
uncertainties shown in Table [21 To the extent that it is 
not possible to judge which model is more physically cor- 
rect, both measurements ought to be treated as equally 
valid. In that situation, the uncertainties, due entirely 
to model assumptions, are roughly ~ 0.4 mag in surface 
brightness and ~ 10% in size. 

7.3. Arp 147 

The #,ST/F814W image of the field Arp 147 contains 
two ring galaxies (Figure [20l Table [3J , one of which has 
a bulge-like component with a tidally disturbed outer re- 
gion (Galaxy 1), and the other is a pure ring (Galaxy 2). 
The best-fitting model for Galaxy 2 is a single-component 
ring, modified by Fourier modes, as seen in Figure l2"0b , 
whereas Galaxy 1 requires two ring components, a bulge, 
and an inner fine-structure component (Figures \20b -h). 
The fine-structure component of Galaxy 1 can be easily 
seen in the surface brightness profile as an upturn within 
r = 0'.'2 of Figure I2TI (left). In addition, the tidal com- 
ponent is slightly bent, which is modeled elegantly using 
the bending modes of Equation [27l As in the case of 
a dust lane, the ring model comes about by truncating 
the inner region of a pure Sersic profile (see Section [5j 
The only difference here is that the truncation radii are a 
larger fraction of the galaxy size. Whereas for the edge- 
on galaxy, it makes more sense to normalize the flux at 
the effective radius (Equation [32]), for ring galaxies, nor- 
malizing the flux at the break radius (Equation 133 j) is 
more intuitive, because it is closer to the peak of the 
profile model. In fact, the peak of the ring is about half- 
way between rb rea k and rb rea k + Ar so f t , but the exact 
location depends on the profile type. 

It is again instructive to compare a traditional fit us- 
ing simple Sersic ellipsoid models (Table [3j bottom) with 
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FlG. 19. — Detailed analysis of an edge-on disk galaxy from GEMS, (a) Original data. (6) Best two-component Sersic profile fit each with 
Fourier modes, corresponding to the parameters shown in Table[2] (c) Best-fit residuals, (d) The fit residuals using traditional (i.e. purely 
ellipsoid) models without masking the dust lane, (e) Residuals after subtracting the best traditional models, masking out the dust lane. 
(J) The bulge component of the best-fit model, (g) The edge-on disk component of the best-fit model, (h) The extended halo component 
of the best-fit model. («) 1-D surface brightness profile. The individual components are shown as dashed lines, and the solid line coursing 
through the data is the sum of the different components. The lower panel shows the residuals of data — model. 
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Table 2. GEMS Disk Galaxy Fitting Results 





# 


— sersic2 — 


Ax ["] 


Ay ["1 


mag/arcsec 2 




r. ["] 


n 


9 


PA [deg] 


Comments 






fouricr 




mode: ampl. & phase [deg] 


mode: ampl. Sz phase [deg] 


mode: ampl. Sz phase [deg] 






# 


— radial — 


Ax ["] 


Ay ["] 




^brcak [ ] 


^ ' soft 


Q 


p PA L UL, &J 




Best 


1 


— sersic2 / — 


0.00 


0.00 


19.80 




0.40 


1.60 


0.72 


44.00 


Trunc. by comp. 


fit 






0.00 


0.00 


0.01 




0.00 


0.01 


0.00 


0.16 


inner: 4 






fouricr 




1: -0.04 


1: -92.16 


3: 


-0.00 


3: -51.36 


4: 0.03 


4: -6.61 












1: 0.00 


1: 3.23 


3 


0.00 


3: 30.65 


4: 0.00 


4: 0.53 








fourier 




5: -0.01 


5: -5.74 


6 


0.01 


6: -7.31 
















5: 0.00 


5: 2.35 


6 


0.00 


6: 1.74 










2 


— sersic2 / — 


{0.00} 


{0.00} 


22.19 




2.29 


0.85 


0.31 


41.30 


Trunc. by comp. 








{0.00} 


{0.00} 


0.02 




0.01 


0.01 


0.00 


0.03 


inner: 4 






fouricr 




1: -0.04 


1: -24.97 


3 


0.02 


3: 26.37 


4: -0.02 


4: -0.04 












1: 0.00 


1: 1.82 


3 


0.00 


3: 1.03 


4: 0.00 


4: 0.79 








fouricr 




5: 0.00 


5: 6.19 


6: 


-0.01 


6: 6.62 
















5: 0.00 


5: 5.63 


6 


0.00 


6: 1.35 










3 


— sersic2 / — 


{0.00} 


{0.00} 


23.92 




4.45 


1.08 


0.49 


41.06 


Trunc. by comp. 








{0.00} 


{0.00} 


0.03 




0.05 


0.02 


0.00 


0.10 


inner: 4 






fouricr 




1: 0.04 


1: 4.59 


3 


0.01 


3: -3.37 


4: -0.01 


4: 40.55 












1: 0.00 


1: 1.16 


3 


0.00 


3: 1.99 


4: 0.00 


4: 3.94 








fouricr 




5: 0.01 


5: -7.19 


6: 


-0.01 


6: -3.25 
















5: 0.00 


5: 1.40 


6 


0.00 


6: 0.73 










4 


— radial — 


0.02 


-0.14 






1.48 


1.48 


0.09 


41.38 


Truncates comp. 








0.00 


0.00 






0.01 


0.02 


0.00 


0.05 


inner: 12 3 






fouricr 




1: -0.12 


1: -156.83 


3 


0.10 


3: -20.62 


4: 0.19 


4: 9.20 












1: 0.00 


1: 1.23 


3 


0.00 


3: 0.52 


4: 0.00 


4: 0.15 








fouricr 




5: 0.10 


5: 2.17 


6 


0.14 


6: 6.62 
















5: 0.00 


5: 0.19 


6 


0.00 


6: 0.10 












merit 


X 2 = 1474348.38 


JV dof = 1435846 




AT free = 64 




= 1.03 




Tradit. 


1 


— scrsic2 — 


0.00 


0.00 


19.98 




0.38 


1.32 


0.74 


42.75 




ellipsoid 






0.00 


0.00 


0.00 




0.00 


0.01 


0.00 


0.31 




model 


2 


— sersic2 — 


{0.00} 


{0.00} 


22.78 




2.57 


0.85 


0.25 


41.21 




with 






{0.00} 


{0.00} 


0.03 




0.02 


0.01 


0.00 


0.06 




dust 


3 


— sersic2 — 


{0.00} 


{0.00} 


23.12 




3.34 


1.77 


0.49 


41.30 




masking 






{0.00} 


{0.00} 


0.04 




0.04 


0.03 


0.00 


0.08 








merit 


X 2 = 1478767.62 


AT dof = 1434511 




AT free = 18 


xl 


= 1.03 





Note. — Best-fitting parameters for an edge-on disk galaxy in GEMS. Sec Table for details. The curly braces ({...}) around parameters 
indicate that they are coupled relative to the first component. Note that the flux amplitude of sersic2 is normalized to the surface brigtness at r e , 
as defined in Equation |32| The "Best fit" parameters (top section) correspond to Panel (b) in Figure fl9l "Traditional ellipsoid model" parameters 
(bottom section) produce residuals shown in Panel (c), and the model is not shown. 
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Fig. 20. — Detailed analysis of Arp 147. (a) Original data. (6) Best Scrsic profile fits of the two galaxies, all with Fourier modes, 
corresponding to the parameters shown in Table \3\ (c) Best-fit residuals, (d) The fit residuals using traditional, i.e. axisymmetric 
ellipsoidal model components, (e) The bulge component of the right-hand galaxy in Panel (6). (J) The inner fine-structure component of 
the best-fit model. (<?) The ring component of the best-fit model, (h) The extended tidal-feature-like component of the best-fit model, (i) 
1-D surface brightness profile of the two galaxies. The individual components are shown as dashed lines, and the solid line coursing through 
the data is the sum of the different components. The lower panel shows the residuals of data — model. 
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Table 3. Arp 147 Fitting Results 





# 


— sersic — 


Ax ["] 


Ay ["] 


mag 


re ["] 


n 


9 


Spa [deg] 


Comments 




# 


— sersic3 — 


Ax ["] 


Ay ["] 


mag/arcsec 2 


re ["] 


n 


q 


P A deg 








fouricr 





mode: ampl. 


& phase [deg] 


mode: ampl. 


& phase [deg] 


mode: ampl. & phase [dcg] 








bending 







mode: 


amplitude ["] 


mode: 


amplitude ["] 


mode: amplitude ["] 






# 


— radial — 












^"forGajk [ ] 


Ar so f t 


'/ 


fPA [ClCgJ 




Galaxy 1 


1 


— sersic — 


0.00 




0.00 


15.49 


1.26 


0.47 


0.50 


194.70 


Bulge. 








0.00 




0.00 


0.00 


0.00 


0.00 


0.00 


0.05 


Inner fine 




2 


— sersic — 


0.07 




0.13 


17.39 


0.35 


0.43 


0.60 


150.75 


structure. 








0.00 




0.00 


0.01 


0.00 


0.01 


0.00 


0.39 






3 


— sersic3 / — 


-0.08 




0.33 


23.84 


1.21 


0.96 


0.18 


187.42 


Trunc. by comp. 






0.01 




0.03 


0.02 


0.02 


0.01 


0.00 


0.03 


inner: 5 (ring). 






fouricr 




1: 


0.03 


1: 49.11 


3: 0.02 


3: 15.69 


4: 0.04 


4: 4.24 










1: 


0.00 


1: 9.08 


3: 0.00 


o. 1.U4 


4: 0.00 


4: 0.45 






4 


— sersic3 / — 


0.09 




-0.34 


22.04 


3.81 


1.94 


0.42 


184.08 


Trunc. by comp. 








0.01 




0.02 


0.01 


0.04 


0.03 


0.00 


0.04 


inner: 5 (Tidal 






fouricr 




1: 


0.16 


1: 21.49 


3: 0.09 


3: 18.36 


4: -0.01 


4: 16.13 


feature) . 










1: 


0.00 


1: 0.82 


3: 0.00 


3: 0.17 


4: 0.00 


4: 1.55 








bending 




2: - 


-0.14 
















5 


— radial — 




2: 


0.00 




10.94 


6.00 


0.18 


187.61 


Truncates comp. 
















0.02 


0.06 


0.00 


0.03 


inner: 3 4 






fouricr 




1: 


0.05 


1: 39.00 


3: 0.02 


3: 19.88 


4: 0.04 


4: 3.88 












1: 


0.00 


1: 4.97 


3: 0.00 


3: 1.18 


4: 0.00 


4: 0.39 




Galaxy 2 


6 


— sersic3 / — 


-18.29 




-7.93 


22.14 


0.78 


1.85 


0.79 


187.91 


Trunc. by comp. 






0.01 




0.01 


0.00 


0.01 


0.01 


0.00 


0.12 


inner: 7 (Ring) 






fouricr 




1: 


0.23 


1: -113.97 


3: 0.07 


3: 15.27 


4: -0.02 


4: 23.33 


magtot = 14.90 










1: 


0.00 


1: 0.53 


3: 0.00 


3: 0.12 


4: 0.00 


4: 0.39 






7 


— radial — 










10.77 


6.08 


0.82 


195.16 


Truncates comp. 
















0.01 


0.01 


0.00 


0.19 


inner: 6 






fourier 
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Note. — Best-fitting parameters for Arp 147. See Table for details. Note that the flux amplitude of sersic3 is normalized to the surface 
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more sophisticated analysis (Table |3l top). In terms of 
the total flux for Galaxy 1, the magnitude of the most so- 
phisticated model is m — 14.18, compared to m — 14.32 
for a model based on classical ellipsoids. Interestingly, 
a single-component fit (not shown) to Galaxy 1 yields 
a magnitude of m = 14.15. For Galaxy 2, we know 
from the outset that classical ellipsoid models are en- 
tirely inappropriate to use. Yet, despite every reason 
to believe that the photometry would be inaccurate, we 
find that the total flux of the traditional ellipsoid fit 
is only 0.2 magnitude different from the most realis- 
tic ring model. These two examples show once again 
that a single-component Sersic ellipsoid fit to compli- 
cated galaxies can produce quite accurate measurement 
of the total flux. 

It is sometimes desirable to conduct bulge-to-disk 
(B/D) decompositions, and Galaxy 1 is an ideal candi- 
date to conduct a comparison. In the traditional ellipsoid 
model (Table El bottom), the B/D ratio is 0.65. The 
more sophisticated model (Table El top) requires sum- 
ming the ring+tidal feature components to obtain the 
disk component, which yields 14.65 mag, thus a B/D ra- 
tio of 0.54. In this situation, most of the differences arise 
from measuring the disk component, which differs by 0.2 
mag, whereas the bulge component is quite robust, with 
a difference of only 0.01 mag. 

It is also of interest to understand how the structural 
parameters are affected by different model choices, in 
particular for the ring Galaxy 2. Whereas the effective 
radius for the ring model is only 0'.'78, for the ellipsoid 
model it is 8'.' 28. This is understandable, bearing in mind 
that the ring has a radius of nearly 8". To a classical 
Sersic profile, the galaxy appears to have a very fiat (in 
fact, a deficit) core, which leads to a low Sersic index of 
n = 0.12. As most of the flux is at 8", beyond which the 
ring flux quickly fades, the ring radius is closely related 
to the effective radius for a classical Sersic model. For 
the inner-truncated ring model (component 7), however, 
the physical size of the ring is captured by the break ra- 
dius rbreak parameter, whereas the r e term no longer has 
the classical meaning of the effective radius (i.e. half the 
light is within r e ). Instead, r e for component 7 is essen- 
tially an exponentially declining scale length parameter, 
given by Equation [33] As the flux dies away quickly be- 
yond the peak, as shown in Figure l2"0"lt. the scale length 
r e for the ring model must therefore be quite small. The 
differences in the r e parameter between the traditional 
model and the truncated model are therefore only due to 
definitions, and not due to systematic or random mea- 
surement uncertainties. 

7.4. M51 

The classical Whirlpool galaxy is a beautiful system 
where a grand-design spiral, M51A, is interacting with 
another spiral, M51B (Figure ED Table B}. In addition 
to there being obvious spiral structures for both galax- 
ies, there are large tidal disturbances that emerge from 
M51B, as seen in the SDSS r-band image provided by 
D. Finkbeiner. Because they are closely overlapping, a 
desirable goal is to deblend M51A and B, as well as to 
model the spiral and tidal structures, simultaneously. 

As with previous examples, we fit this galaxy using 
both the most sophisticated analysis (Table HI top) in 
our toolbox, and comparing the results to the traditional 



axisymmetric ellipsoids (Table |H bottom) analysis. The 
traditional analysis requires two components each, in or- 
der to decompose a galaxy ostensibly into a bulge and a 
disk. The reduction in xt between the two methods is 
modest, because most of the residuals come from high- 
frequency starforming regions that are not removed by 
models which are fundamentally smooth, despite being 
modified by radial Fourier modes and spiral rotations. 

In the most detailed analysis of M51A, we use two spi- 
ral arm components and two components for the bulge. 
There is actually not a strong need to use two compo- 
nents for the bulge except to better capture the detailed 
profile shape, which has an inflection at r f=s 0'A, as seen 
in Figure |2"TK On the other hand, the use of two spi- 
ral components is necessary because the spiral arm has 
a "kink" in the rotation that cannot be created by us- 
ing a single smooth rotation function. The spiral struc- 
tures are modified by Fourier modes to create both a 
slight lopsidedness and other subtle features. Because 
there are more degrees of freedom in a two-arm spiral, 
the higher order Fourier modes also can "see" detailed 
structures, like the reverse flaring of the spiral structure 
in Figure [2T]f . 

For M51B, we employ three components in the fit, a 
bulge (component 5 in Tabled top), a tidal feature com- 
ponent (component 6), and a spiral function (component 
7), which model the three most visually striking compo- 
nents. The tidal feature is mostly obtained by using the 
second and third bending modes of Equation [27l as il- 
lustrated in Figure 1101 However, bending modes 2 and 
3 are symmetric functions, so the high degree of asym- 
metry comes about because of combined action with the 
Fourier modes, which is shown to have a high amplitude 
of 0.23 for the m = 1 mode, as well as moderate values 
for other modes. Incidentally, despite the complexity of 
the higher order structures, all the parameter values are 
determined automatically by Galfit without the need 
for an user to provide initial guesses (i.e. initially all 
values) and without tweaking at any point in the analysis 
(which is hardly feasible anyhow). 

For even those who are experienced with detailed para- 
metric fitting, one of the alarming facts about this anal- 
ysis is that it employs 103 free parameters in the best- 
fit model. So there are natural concerns about param- 
eter degeneracies. However, as we have discussed in 
Section 16.51 parameter degeneracies do not arise purely 
based on the number of free parameters, but rather on 
the types of parameters involved. The availability of 
spatial information in 2-D provides one of the most im- 
portant ways to break parameter degeneracies. We see 
this explicitly in Figures l2"l"ls-/t. where there is little evi- 
dence that the subcomponents for M51A are strongly in- 
fluenced by M51B, and vice versa. Furthermore, within 
each galaxy, the subcomponents are so different in shape, 
both qualitatively and quantitatively, that the amount of 
crosstalk between them is also not significant. Therefore, 
despite the extreme complexity of this system, and the 
use of 103 free parameters, we find that degeneracies be- 
tween the parameters are not an issue. Or, if they exist, 
they do so at a low enough level that they do not sig- 
nificantly affect the main parameters of interest, like the 
luminosity of the subcomponents, or the profile shapes 
and sizes. 

There are, however, seemingly degenerate conditions 



34 




tfi © If, 

( ii i in .1.1 n i v 



[ 3 oasojB / §buu gv] Jssas ri 



■tf C\2 OCAi 

do do 
i i 



Fig. 21. — Detailed analysis of M51. (a) Original data. (6) Best Sersic profile fits of the M51A and B, all with Fourier modes, 
corresponding to the parameters shown in Table [4] (c) Best- fit residuals, (d) The fit residuals using traditional, axisymmetric, ellipsoidal 
model components, (e) The extended grand-design spiral component of M51A model in Panel (6). (/) The inner fine-structure spiral 
component of the best-fit model, (g) The spiral component of M51B. (h) The extended tidal feature-like component of M51B, using 
simultaneous bending and Fourier modes. A bulge component is present but not shown in the figures of M51A and B. i) 1-D surface 
brightness profile of the two galaxies. The individual components are shown as dashed lines, and the solid line coursing through the data 
is the sum of the different components. The lower panel shows the residuals of data — model. 
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Table 4. M51 Fitting Results 
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that have little to do with parameter coupling. Instead, 
these are attributed to the fact that M51 has many of 
non-smooth structural features, caused by dust lanes, 
star-forming regions, tidal disturbances, and so forth. 
Such spatially localized features, if strong enough, can 
influence Galfit to "lock" on to them if the initial 
conditions happened to be sufficiently close. The con- 
sequences appear as degeneracies when, in fact, there 
are many small local minima solutions. This graininess 
in the % 2 terrain introduces slight perturbations to the 
models, and may even cause fairly large shape differ- 
ences in the final solutions. However, to a large extent, 
it rarely affects the main parameters of interest, such 
as the luminosity of a particular component or its size, 
which are determined by much more global features than 
the nuisances of local fluctuations to which higher order 
parameters are more sensitive. 

To gain some intuitive insight into the effects of com- 
plex analysis, it is instructive to compare simple and 
complex methods with regard to global and subcompo- 
nent properties. In terms of the total luminosity, here we 
find excellent agreement between sophisticated and tradi- 
tional analysis, respectively, of m r — 8.24 vs. m r = 8.25 
for M51A, and m r = 8.80 vs. m r = 8.73 for M51B. While 
this level of agreement may at first seem surprising, it is 
expected given the basic premise of least-squares mini- 
mization. In fact, even a single-component fit to M51A 
yields m r — 8.0, and for M51B m r — 9.0, which are 
both quite close to the overall best-fit models, despite 
the complications in the image. The main reason for 
the discrepancy here is the uncertainty in the sky, due 
to there being a large gradient. This fundamentally sets 
the limit on the accuracy of the photometry to perhaps 
no better than 0.1 to 0.2 magnitude, independent of the 
analysis method. 

The most sensitive benchmark for understanding dif- 
ferences in the analysis is in detailed decompositions. 
Here we compare the bulge-to-disk decomposition re- 
sults. In the traditional ellipsoid analysis, we find a B/D 
ratio of 0.23 for M51A and 4.9 for M51B. The large B/D 
ratio for M51B is clearly unphysical, and is driven by 
the large Sersic index (n = 8.0) of the bulge component, 
which is increased to accommodate the flux in the out- 
skirts due to tidal features. In the most detailed analysis, 
the B/D ratio for M51A is 0.16, whereas for M51B it is 
merely 0.17. Examining the bulge of M51A more closely, 
we find that the detailed analysis yields a total flux of 
10.38 mag, whereas the traditional analysis extracts a 
brighter bulge of 10.05 mag. The differences come from 
the fact that the light of the inner spiral is in part driving 
up the Sersic index of the bulge when it is not properly 
accounted. It is probably safe to conclude that a magni- 
tude of 10.05 is a firm upper limit to the bulge luminosity. 

Finally, it is worthwhile to compare how the disk pa- 
rameters differ between the analyses to gain an under- 
standing for how coordinate rotation affects the inter- 
pretation of the parameters for the spiral models. From 
Table |4j we find that the Sersic index of the simple and 
complex models are essentially identical for M51A, at 
n 0.33. The interpretation for M51B is more compli- 
cated, because the "disk" in an ellipsoidal model is not 
qualitatively the same structure as the spiral analysis. In 
fact, it is necessary to hold the Sersic index of the tidal 
component 6 fixed in the analysis. Nevertheless there are 



clearly quantitative differences in that the simple anal- 
ysis is larger by 55% in n. With regard to the effective 
radius, the traditional analysis of M51A finds the disk 
size to be about 2.'2, which compares favorably with the 
spiral model size of 2.'8, or a 25% difference. Further- 
more, the disk magnitudes for M51A differ only by 0.03 
mag between simple and complex. 

These comparisons therefore demonstrate that despite 
the complex analysis being much more realistic looking, 
fundamentally the meaning of the structural parameters 
(size, luminosity, concentration index) are unchanged 
from the original definition, even in the situation of spi- 
ral components. This is an useful fact because our prior 
intuitions, honed on fitting ellipsoidal models, continue 
to be applicable. We note that the generally good agree- 
ment between detailed and simplistic analysis witnessed 
here and in previous examples is not entirely coinciden- 
tal. It so happens because all shapes are fundamentally 
perturbations of the best-fitting ellipsoidal model, even 
if the result bears no resemblance to the original ellipse. 

7.5. NGC 289 

NGC 289 is an SAB(rs)bc galaxy, with a weak bar 
and a complex inner spiral system (Figure Table [SJ 
that resembles a ring. Upon closer examination, the ring 
appearance comes about because there exists a bifurca- 
tion in the spiral structure that connects up with the 
opposing spiral arm. Furthermore, the bar is also multi- 
component, with the inner component oriented at an an- 
gle nearly 45° from the strong outer bar. 

The best-fit analysis involves three spiral components, 
an inner and an outer bar component, and a bulge (Ta- 
ble[5j top). All except for the bulge component are modi- 
fied by five Fourier modes, and are shown in Figures I2"2"1s- 
h. The requirement of components 3 and 4 (Figures |2"2T 
and g) is clear, because they are what form the most 
striking and intricate patterns in the center, while the 
requirement of component 5 (Figure 122f t) is only evident 
in the residuals, and makes up some of the diffuse light 
within the inner 60" region. Although it does not seem 
like an essential component, the inner bar structure (Fig- 
urel2"2"le) qualitatively affects the detailed residual pattern 
at the center, and is therefore included. When all the de- 
tailed inner structures are properly accounted for, it is 
straightforward to infer the bulge component, and assess 
the uncertainties by varying different parameters of the 
bulge. Doing so does not affect the inner fine-structures 
because they are sharp and well localized. 

Conducting the same decomposition using traditional 
ellipsoid models (Table [5j bottom), we opted to fit three 
components, ostensibly to model a bulge, disk, and a bar. 
The result produces residuals seen in Figure |2"2"H . reveal- 
ing the intricate details of the inner spiral system. From 
the fit, even though the disk and the bar component are 
sensible, the bulge component is actually fitting a dif- 
fuse disk component, which, in retrospect, is that shown 
in Figure \22h . Because that inner spiral component is 
quite luminous, and because there exists a bulge com- 
ponent superposed on top of it, this quasi-bulge model 
is almost 0.7 mag brighter than that inferred through 
the detailed modeling above. Adding a fourth ellipsoid 
model is not possible, because the central spiral residuals 
are so great that they completely suppress the addition 
of another component, causing the flux to go to zero. 
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Fig. 22. — Detailed analysis of NGC 289 from CINGS. (a) Original data. (6) Best Sersic profile fits with spiral rotation functions 
and Fourier modes, corresponding to the parameters shown in Table (c) Best-fit residuals, (d) The fit residuals using traditional, 
axisymmctric, ellipsoidal model components, (e) The fine details of the inner bar structure of Panel (b). (J) Spiral component 1 of 3 of 
the best-fit model, (g) The spiral component 2 of 3. (h) The spiral component 3 of 3. A bulge component is present but not shown in the 
figures, (i) 1-D surface brightness profile of the galaxy. The individual components are shown as dashed lines, and the solid line coursing 
through the data is the sum of the different components. The lower panel shows the residuals of data — model. 



Table 5. NGC 289 Fitting Results 
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Once again, comparing the total luminosity between 
the best-fit model with the ellipsoid fit, we find an ex- 
cellent agreement of m = 10.37 mag vs. m = 10.42, 
respectively. For a single-component ellipsoid fit, there 
is also an excellent agreement of m = 10.46, despite the 
main structural details not being unaccounted. 

8. DISCUSSION AND CONCLUSIONS 

This study is a proof of concept for how to conduct 
more realistic image-fitting analysis using purely para- 
metric functions, by breaking free from traditional as- 
sumptions about axisymmctry. We introduced several 
new ideas, including the use of Fourier and bending 
modes, spiral rotation functions, and truncation func- 
tions. These features can be used individually, or com- 
bined in arbitrary ways. While these features are individ- 
ually quite simplistic, used collectively they proliferate a 
dizzying array of possibilities. Even so, the interpreta- 
tion of each component remains intuitive, down to the 
meaning of each fitting parameter. Indeed, the interpre- 
tation of the traditional ellipsoidal profile parameters, 
such as those for the Sersic function, remains essentially 
unchanged under modification. We then applied these 
techniques to five case studies, illustrating that highly 
complex and intricate structures can be modeled using 
fully parametric techniques. 

There are many practical applications for these tech- 
niques. For instance, the Fourier modes are useful for 
quantifying the average global symmetry of a galaxy, 
and can easily be automated for galaxy surveys. It is 
also possible to disentangle bright from faint asymme- 
tries, and to conduct more robust bulge-to-disk decom- 
positions in some galaxies. It would be useful to quantify 
how much of the total flux is in a bulge versus the tidally 
distorted component, which has implications for issues 
such as late- vs. early-stage mergers, or major vs. minor 
mergers. 

More than just a presentation of new techniques, one of 
the main purposes of this study is to highlight a method 
to more realistically quantify measurement uncertainties 
in high-S/N images. In galaxy fitting, the most desirable 
goal will always be to obtain a fit with the lowest x 2 > us- 
ing the simplest model. In the past, this idea was closely 
tied to the practice of using one- or two-component el- 
lipsoid models, out of necessity. Simplicity is not nec- 
essarily congruent with propriety or reality. This study 
promotes the notion that simple models can be realistic. 
It also opens up new possibilities for more detailed anal- 
ysis depending on the image complexity. However, this 
possibility is both a blessing and a curse. For, the fact 
there is not one generic solution for any galaxy leads to 
the following conundrum in image analysis, but one that 
illustrates the merit of our approach: 

"What model should one adopt, how much detail is 
enough, and what about degeneracies?" We have 

shown that detailed decomposition analysis can be arbi- 
trarily sophisticated. It is for that same reason there is 
often not a single, unique answer. However, the essential 
fact, as seen through our examples and other detailed 
analysis outside of this work, is that the marginal return 
of adding complexity quickly diminishes. Therefore, the 
above conundrum is in practice easy to address by con- 
ducting analyses of varying sophistication without prej- 
udice, then judging the outcome by taking a clear view 



of what goal is to be achieved. If different solutions yield 
the same result for a desired science goal (e.g., bulge 
luminosity, bulge-to-disk ratio, average size, total lumi- 
nosity, etc.), then it does not matter which model one 
adopts. If they yield different outcomes, then the most 
realistic analysis ought to be the more true. However, if 
it is not possible to decide on the correct model, the dif- 
ferent results by definition give an estimate of the model- 
dependent measurement uncertainty. This last attribute, 
rather than being a perceived weakness, is fundamentally 
that which makes the analysis quantitatively rigorous. 

Despite the flexibility allowed by the models, this pa- 
per is merely an initial demonstration of concept and 
leaves many issues unsolved. Currently, the formulation 
of the spiral rotation function is fairly rigid, and can- 
not produce arms that wind back onto itself (although 
that can be approximated using the ring feature in Gal- 
fit) . The amount of curvature in the bending modes can 
only fit arcs and not fuller semi-circles (which can partly 
be modeled using a lopsided ring). There remains sub- 
stantial room for future growth in profile types, shape 
definitions, and toward spatial-spectral decompositions 
for integral-field data. 
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Fig. 23. — Appendix Figure A: Schematics of a hyperbolic tangent rotation function. n n is the inner radius where the rotation angle 
reaches 20 degrees relative to the PA of the best fitting ellipse of a component. Below r- ln the function flattens out to zero degrees. r ou t is 
the outer radius, beyond which the function flattens out to a constant rotation angle, and 8 ou t is the total amount of rotation out to r Q ut- 
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APPENDIX 

HYPERBOLIC TANGENT ROTATION FUNCTION 



The hyperbolic tangent (tanh(r 



1 Tout ) 



(9pA' r )) portion of the a-tanh (Equation |28|) and log-tanh (Equation |29|) 



rotation functions is given by Equation [5] below. The constant CDEF is defined such that at the mathematical "bar 
radius" r- m , the rotation angle 9 reaches 20°. This definition is entirely empirical. The above Figure shows a pure tanh 
rotation function, where the rotation angle reaches a maximum 9 out near r = r out . Beyond r out , the rotation angle 
levels off at f? out . This function is multiplied with either a logarithmic or a power-law function to produce the desired 
asymptotic rotation behavior seen in more realistic galaxies (see Section [J). 



CDEF = 0.23 (constant for "bar" definition) 



(1) 



_ 2 x CDEF 
~ |0 O ut| + CDEF 

B = (2- tanh" 1 (A)) 



- 1.00001 

Tout 
'"out 'in 



= yj Ax 2 + Ay 2 (circular — centric distance) 



tanh(r in! r ou t,^inci,^PA'; r ) = 0.5 x ( tanh 



s (^-i 

\ r out 



(2) 

(3) 
(4) 

(5) 
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B — HYPERBOLIC TANGENT TRUNCATION FUNCTION 



The hyperbolic tangent truncation function (tanh(a;o, yo; ?"break; ^soft> 1, #pa)) (see Section [5]) is very similar to the 
coordinate rotation formulation in Appendix A, except for different constants that define the flux ratio at the truncation 
radii: at r = rb rea k the flux is 99% of the untruncated model profile, whereas at r = r so f t the flux is 1%. With this 
definition, Equation [7] is the truncation function: 



B = 2.65-4.98 



''"break 



\ ^brcak ^*soft 



tanh(x , yo] r-brcak, r soft , q, 9 PA ) = 0.5 x tanh 



(2-B) 



B 



'break 



i 



(6) 



Note that the radius r is a generalized radius (as opposed to a circular-centric distance), i.e. one that is perturbed by 
Co, bending modes, or Fourier modes, of the truncation function. When the softening length (Ar so f t ) is used as a free 
parameter, it is defined as Ar so f t = rbreak — ?"soft- 
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